From 165468fcceaf7d52c2fd24cbf668bfd94564ae9f Mon Sep 17 00:00:00 2001
From: Arne Naegel <naegel@gcsc.uni-frankfurt.de>
Date: Fri, 14 Jun 2024 07:46:49 +0200
Subject: [PATCH 1/4] Protecting output

---
 content/brom/BromDiffusion.ipynb | 262 +++----------------------------
 content/brom/BromDiffusion.py    |  30 ++--
 2 files changed, 34 insertions(+), 258 deletions(-)

diff --git a/content/brom/BromDiffusion.ipynb b/content/brom/BromDiffusion.ipynb
index 2c2753e..3bac8f2 100644
--- a/content/brom/BromDiffusion.ipynb
+++ b/content/brom/BromDiffusion.ipynb
@@ -20,7 +20,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "e6150636",
    "metadata": {},
    "outputs": [],
@@ -43,22 +43,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "id": "42847742",
    "metadata": {},
-   "outputs": [
-    {
-     "ename": "AttributeError",
-     "evalue": "module 'ug4py.pyugcore' has no attribute 'Domain3d'",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[3], line 9\u001b[0m\n\u001b[1;32m      5\u001b[0m numRefs \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m  \u001b[38;5;66;03m# Number of Refinement steps on given grid\u001b[39;00m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;66;03m# Choosing a domain object\u001b[39;00m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m# (either 1d, 2d or 3d suffix)\u001b[39;00m\n\u001b[0;32m----> 9\u001b[0m dom \u001b[38;5;241m=\u001b[39m \u001b[43mug4\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mDomain3d\u001b[49m()\n\u001b[1;32m     11\u001b[0m \u001b[38;5;66;03m# Loading the given grid into the domain\u001b[39;00m\n\u001b[1;32m     12\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLoading Domain \u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m+\u001b[39m gridName \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m...\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
-      "\u001b[0;31mAttributeError\u001b[0m: module 'ug4py.pyugcore' has no attribute 'Domain3d'"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Setup:\n",
     "# defining needed subsets, grid and number of refinements\n",
@@ -86,7 +74,7 @@
     "\n",
     "    print(\"Refining done\")\n",
     "\n",
-    "# checking if geometry has the needed subsets of the probelm\n",
+    "# checking if geometry has the needed subsets of the problem\n",
     "sh = dom.subset_handler()\n",
     "for e in requiredSubsets:\n",
     "    if sh.get_subset_index(e) == -1:\n",
@@ -105,22 +93,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "id": "8f75da00",
    "metadata": {},
-   "outputs": [
-    {
-     "ename": "AttributeError",
-     "evalue": "module 'ug4py.pyugcore' has no attribute 'ApproximationSpace3d'",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[4], line 6\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;28mtype\u001b[39m \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mLagrange\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m      4\u001b[0m order \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m  \u001b[38;5;66;03m# polynom order for lagrange\u001b[39;00m\n\u001b[0;32m----> 6\u001b[0m approxSpace \u001b[38;5;241m=\u001b[39m \u001b[43mug4\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mApproximationSpace3d\u001b[49m(dom)\n\u001b[1;32m      7\u001b[0m approxSpace\u001b[38;5;241m.\u001b[39madd_fct(fct, \u001b[38;5;28mtype\u001b[39m, order)\n\u001b[1;32m      8\u001b[0m approxSpace\u001b[38;5;241m.\u001b[39minit_levels()\n",
-      "\u001b[0;31mAttributeError\u001b[0m: module 'ug4py.pyugcore' has no attribute 'ApproximationSpace3d'"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "\n",
     "# Create approximation space which describes the unknowns in the equation\n",
@@ -181,7 +157,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "id": "f8e0894b",
    "metadata": {},
    "outputs": [],
@@ -213,7 +189,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "id": "3140f477",
    "metadata": {},
    "outputs": [],
@@ -255,179 +231,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "id": "1b79efbe",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "+++ Integrating: [0, 1]\n",
-      "+++ Timestep +++1\n",
-      "+++ Reassemble (t=0, dt=0.00625)\n",
-      "+++ Timestep +++2\n",
-      "+++ Timestep +++3\n",
-      "+++ Timestep +++4\n",
-      "+++ Timestep +++5\n",
-      "+++ Timestep +++6\n",
-      "+++ Timestep +++7\n",
-      "+++ Timestep +++8\n",
-      "+++ Timestep +++9\n",
-      "+++ Timestep +++10\n",
-      "+++ Timestep +++11\n",
-      "+++ Timestep +++12\n",
-      "+++ Timestep +++13\n",
-      "+++ Timestep +++14\n",
-      "+++ Timestep +++15\n",
-      "+++ Timestep +++16\n",
-      "+++ Timestep +++17\n",
-      "+++ Timestep +++18\n",
-      "+++ Timestep +++19\n",
-      "+++ Timestep +++20\n",
-      "+++ Timestep +++21\n",
-      "+++ Timestep +++22\n",
-      "+++ Timestep +++23\n",
-      "+++ Timestep +++24\n",
-      "+++ Timestep +++25\n",
-      "+++ Timestep +++26\n",
-      "+++ Timestep +++27\n",
-      "+++ Timestep +++28\n",
-      "+++ Timestep +++29\n",
-      "+++ Timestep +++30\n",
-      "+++ Timestep +++31\n",
-      "+++ Timestep +++32\n",
-      "+++ Timestep +++33\n",
-      "+++ Timestep +++34\n",
-      "+++ Timestep +++35\n",
-      "+++ Timestep +++36\n",
-      "+++ Timestep +++37\n",
-      "+++ Timestep +++38\n",
-      "+++ Timestep +++39\n",
-      "+++ Timestep +++40\n",
-      "+++ Timestep +++41\n",
-      "+++ Timestep +++42\n",
-      "+++ Timestep +++43\n",
-      "+++ Timestep +++44\n",
-      "+++ Timestep +++45\n",
-      "+++ Timestep +++46\n",
-      "+++ Timestep +++47\n",
-      "+++ Timestep +++48\n",
-      "+++ Timestep +++49\n",
-      "+++ Timestep +++50\n",
-      "+++ Timestep +++51\n",
-      "+++ Timestep +++52\n",
-      "+++ Timestep +++53\n",
-      "+++ Timestep +++54\n",
-      "+++ Timestep +++55\n",
-      "+++ Timestep +++56\n",
-      "+++ Timestep +++57\n",
-      "+++ Timestep +++58\n",
-      "+++ Timestep +++59\n",
-      "+++ Timestep +++60\n",
-      "+++ Timestep +++61\n",
-      "+++ Timestep +++62\n",
-      "+++ Timestep +++63\n",
-      "+++ Timestep +++64\n",
-      "+++ Timestep +++65\n",
-      "+++ Timestep +++66\n",
-      "+++ Timestep +++67\n",
-      "+++ Timestep +++68\n",
-      "+++ Timestep +++69\n",
-      "+++ Timestep +++70\n",
-      "+++ Timestep +++71\n",
-      "+++ Timestep +++72\n",
-      "+++ Timestep +++73\n",
-      "+++ Timestep +++74\n",
-      "+++ Timestep +++75\n",
-      "+++ Timestep +++76\n",
-      "+++ Timestep +++77\n",
-      "+++ Timestep +++78\n",
-      "+++ Timestep +++79\n",
-      "+++ Timestep +++80\n",
-      "+++ Timestep +++81\n",
-      "+++ Timestep +++82\n",
-      "+++ Timestep +++83\n",
-      "+++ Timestep +++84\n",
-      "+++ Timestep +++85\n",
-      "+++ Timestep +++86\n",
-      "+++ Timestep +++87\n",
-      "+++ Timestep +++88\n",
-      "+++ Timestep +++89\n",
-      "+++ Timestep +++90\n",
-      "+++ Timestep +++91\n",
-      "+++ Timestep +++92\n",
-      "+++ Timestep +++93\n",
-      "+++ Timestep +++94\n",
-      "+++ Timestep +++95\n",
-      "+++ Timestep +++96\n",
-      "+++ Timestep +++97\n",
-      "+++ Timestep +++98\n",
-      "+++ Timestep +++99\n",
-      "+++ Timestep +++100\n",
-      "+++ Timestep +++101\n",
-      "+++ Timestep +++102\n",
-      "+++ Timestep +++103\n",
-      "+++ Timestep +++104\n",
-      "+++ Timestep +++105\n",
-      "+++ Timestep +++106\n",
-      "+++ Timestep +++107\n",
-      "+++ Timestep +++108\n",
-      "+++ Timestep +++109\n",
-      "+++ Timestep +++110\n",
-      "+++ Timestep +++111\n",
-      "+++ Timestep +++112\n",
-      "+++ Timestep +++113\n",
-      "+++ Timestep +++114\n",
-      "+++ Timestep +++115\n",
-      "+++ Timestep +++116\n",
-      "+++ Timestep +++117\n",
-      "+++ Timestep +++118\n",
-      "+++ Timestep +++119\n",
-      "+++ Timestep +++120\n",
-      "+++ Timestep +++121\n",
-      "+++ Timestep +++122\n",
-      "+++ Timestep +++123\n",
-      "+++ Timestep +++124\n",
-      "+++ Timestep +++125\n",
-      "+++ Timestep +++126\n",
-      "+++ Timestep +++127\n",
-      "+++ Timestep +++128\n",
-      "+++ Timestep +++129\n",
-      "+++ Timestep +++130\n",
-      "+++ Timestep +++131\n",
-      "+++ Timestep +++132\n",
-      "+++ Timestep +++133\n",
-      "+++ Timestep +++134\n",
-      "+++ Timestep +++135\n",
-      "+++ Timestep +++136\n",
-      "+++ Timestep +++137\n",
-      "+++ Timestep +++138\n",
-      "+++ Timestep +++139\n",
-      "+++ Timestep +++140\n",
-      "+++ Timestep +++141\n",
-      "+++ Timestep +++142\n",
-      "+++ Timestep +++143\n",
-      "+++ Timestep +++144\n",
-      "+++ Timestep +++145\n",
-      "+++ Timestep +++146\n",
-      "+++ Timestep +++147\n",
-      "+++ Timestep +++148\n",
-      "+++ Timestep +++149\n",
-      "+++ Timestep +++150\n",
-      "+++ Timestep +++151\n",
-      "+++ Timestep +++152\n",
-      "+++ Timestep +++153\n",
-      "+++ Timestep +++154\n",
-      "+++ Timestep +++155\n",
-      "+++ Timestep +++156\n",
-      "+++ Timestep +++157\n",
-      "+++ Timestep +++158\n",
-      "+++ Timestep +++159\n",
-      "+++ Timestep +++160\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "\n",
     "# Define start time, end time and step size\n",
@@ -468,53 +275,22 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "c4a48c79",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      "Pyvista input: \n",
-      "UnstructuredGrid (0x1042a86a0)\n",
-      "  N Cells:    5376\n",
-      "  N Points:   3509\n",
-      "  X Bounds:   -1.000e+00, 1.000e+00\n",
-      "  Y Bounds:   -1.000e+00, 1.000e+00\n",
-      "  Z Bounds:   0.000e+00, 3.500e+00\n",
-      "  N Arrays:   1\n",
-      "\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "4830ded197494a0cb69378fe42e07e22",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "Widget(value=\"<iframe src='http://localhost:55551/index.html?ui=P_0x2a0e72a10_1&reconnect=auto' style='width: …"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Exporting the result to a vtu-file\n",
     "# can be visualized in paraview or with a python extension\n",
     "ug4.WriteGridFunctionToVTK(usol, \"BromFinal\")\n",
     "\n",
-    "# Plotting the result using pyvista\n",
-    "import pyvista\n",
-    "result = pyvista.read('BromFinal.vtu')\n",
-    "print()\n",
-    "print(\"Pyvista input: \")\n",
-    "print(result)\n",
-    "result.plot(scalars=\"u\", show_edges=True, cmap='hot')\n"
+    "try:\n",
+    "    # Plotting the result using pyvista\n",
+    "    import pyvista\n",
+    "    result = pyvista.read('BromFinal.vtu')\n",
+    "    result.plot(scalars=\"u\", show_edges=True, cmap='hot')\n",
+    "except:\n",
+    "     print(\"Plotting failed.\")"
    ]
   },
   {
diff --git a/content/brom/BromDiffusion.py b/content/brom/BromDiffusion.py
index 1f0392a..e8f0e72 100644
--- a/content/brom/BromDiffusion.py
+++ b/content/brom/BromDiffusion.py
@@ -13,7 +13,7 @@
 # 
 # ## Setup
 
-# In[2]:
+# In[ ]:
 
 
 import sys
@@ -26,7 +26,7 @@
 
 # ## Grid and domain
 
-# In[3]:
+# In[ ]:
 
 
 # Setup:
@@ -55,7 +55,7 @@
 
     print("Refining done")
 
-# checking if geometry has the needed subsets of the probelm
+# checking if geometry has the needed subsets of the problem
 sh = dom.subset_handler()
 for e in requiredSubsets:
     if sh.get_subset_index(e) == -1:
@@ -66,7 +66,7 @@
 
 # ## Approximation space 
 
-# In[4]:
+# In[ ]:
 
 
 # Create approximation space which describes the unknowns in the equation
@@ -114,7 +114,7 @@
 domainDisc.add(dirichletBND)
 
 
-# In[4]:
+# In[ ]:
 
 
 # Technical stuff
@@ -136,7 +136,7 @@ def MyFunction_v (u,v):
 
 # ## Solve transient problem
 
-# In[5]:
+# In[ ]:
 
 
 # Use the Approximationspace to
@@ -168,7 +168,7 @@ def MyInitialValue(x, y, z, t, si):
 lsolver=ug4.LUCPU1()
 
 
-# In[6]:
+# In[ ]:
 
 
 # Define start time, end time and step size
@@ -201,20 +201,20 @@ def MyInitialValue(x, y, z, t, si):
 
 # ## Exporting the result 
 
-# In[8]:
+# In[ ]:
 
 
 # Exporting the result to a vtu-file
 # can be visualized in paraview or with a python extension
 ug4.WriteGridFunctionToVTK(usol, "BromFinal")
 
-# Plotting the result using pyvista
-import pyvista
-result = pyvista.read('BromFinal.vtu')
-print()
-print("Pyvista input: ")
-print(result)
-result.plot(scalars="u", show_edges=True, cmap='hot')
+try:
+    # Plotting the result using pyvista
+    import pyvista
+    result = pyvista.read('BromFinal.vtu')
+    result.plot(scalars="u", show_edges=True, cmap='hot')
+except:
+     print("Plotting failed.")
 
 
 # In[ ]:

From a0564e3f78dc3deb6c60180338d9b6ee5966b1c5 Mon Sep 17 00:00:00 2001
From: Arne Naegel <naegel@gcsc.uni-frankfurt.de>
Date: Fri, 14 Jun 2024 07:49:02 +0200
Subject: [PATCH 2/4] Adding test for bromide diffusion

---
 tests/integration/test-solvers.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/tests/integration/test-solvers.py b/tests/integration/test-solvers.py
index 6631f1d..9a57cad 100644
--- a/tests/integration/test-solvers.py
+++ b/tests/integration/test-solvers.py
@@ -65,7 +65,6 @@ def main():
 class TestAll(unittest.TestCase):
     def setUp(self):
         os.environ['PYVISTA_OFF_SCREEN'] = "True"
-        set_test_dir("content/tutorial-solver")
 
     def tearDown(self):
         set_test_dir("../..")
@@ -74,10 +73,16 @@ def tearDown(self):
          # self.assertEqual(test_dir("content/tutorial-solver"), True, "Must return true")  
 
     def test_smoothers(self):
+        set_test_dir("content/tutorial-solver")
         self.assertEqual(test_file("example01-smoothers.py"), 0, "Must return true")
     
     def test_multigrid(self):
+        set_test_dir("content/tutorial-solver")
         self.assertEqual(test_file("example03-multigrid.py"), 0, "Must return true")
+    
+    def test_bromide(self):
+        set_test_dir("content/brom")
+        self.assertEqual(test_file("BromDiffusion.py"), 0, "Must return true")
         
 if __name__ == "__main__":
     # main()

From a2c3762d0a91d38046e23249f584751fd13d89c1 Mon Sep 17 00:00:00 2001
From: Arne Naegel <naegel@gcsc.uni-frankfurt.de>
Date: Thu, 27 Jun 2024 09:02:57 +0200
Subject: [PATCH 3/4] Update for skin. Working version for binding.

---
 content/skin/SkinDiffusion.ipynb            |   10 +-
 content/skin/SkinDiffusion.py               |  288 +++--
 content/skin/SkinDiffusionWithBinding.ipynb |  913 +++-----------
 content/skin/SkinDiffusionWithBinding.py    |  219 ++++
 content/skin/SkinDiffusionWithLagtime.ipynb | 1205 ++++++++++++++++++-
 content/skin/SkinDiffusionWithLagtime.py    |  201 ++--
 6 files changed, 1860 insertions(+), 976 deletions(-)
 create mode 100644 content/skin/SkinDiffusionWithBinding.py

diff --git a/content/skin/SkinDiffusion.ipynb b/content/skin/SkinDiffusion.ipynb
index eca6ea7..9cdb041 100644
--- a/content/skin/SkinDiffusion.ipynb
+++ b/content/skin/SkinDiffusion.ipynb
@@ -27,7 +27,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -38,7 +38,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -76,11 +76,11 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Loading Domain {gridName}...\n",
+      "Loading Domain 'skin2d-aniso.ugx'...\n",
       "Domain loaded.\n",
       "Refining ...\n",
-      "Refining step {i} ...\n",
-      "Refining step {i} ...\n",
+      "Refining step {0} ...\n",
+      "Refining step {1} ...\n",
       "Refining done\n"
      ]
     }
diff --git a/content/skin/SkinDiffusion.py b/content/skin/SkinDiffusion.py
index a903169..6a7dfaa 100644
--- a/content/skin/SkinDiffusion.py
+++ b/content/skin/SkinDiffusion.py
@@ -1,146 +1,260 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# [<img src="../../header.svg">](../index.ipynb)
+# 
+# ---
+
+# # Application: Simulation of Drug Transport across a Virtual Skin Membrane
+# 
+# [Advanced version](SkinDiffusionWithLagtime.ipybnd)
+
+# ### Preliminaries
+
+# In[11]:
+
+
 import ug4py.pyugcore as ug4
 import ug4py.pylimex as limex
 import ug4py.pyconvectiondiffusion as cd
-# import ug4py.pysuperlu as slu
 
 
-# Setup:
-# defining needed subsets, gird and number of refinements
+# In[12]:
+
+
+# Importing
+import sys
+sys.path.append('..')
+import modsimtools as util
+
+
+# ## Solve problem
+
+# In[3]:
+
+
+# Setup: Defining requires subsets, grid and number of refinements
 requiredSubsets = ["LIP", "COR", "BOTTOM_SC", "TOP_SC"] # defining subsets
-gridName = "skin2d-aniso.ugx"  # Grid created with ProMesh
+gridName = 'skin2d-aniso.ugx'  # Grid created with ProMesh
 numRefs = 2  # Number of Refinement steps on given grid
 
-# Choosing a domain object
-# (either 1d, 2d or 3d suffix)
-dom = ug4.Domain2d()
-
-# Loading the given grid into the domain
-print(f"Loading Domain {gridName} ...")
-ug4.LoadDomain(dom, gridName)
-print("Domain loaded.")
-
-# Optional: Refining the grid
-if numRefs > 0:
-    print("Refining ...")
-    refiner = ug4.GlobalDomainRefiner(dom)
-    for i in range(numRefs):
-        ug4.TerminateAbortedRun()
-        refiner.refine()
-        print(f"Refining step {i} ...")
-
-    print("Refining done")
-
-# checking if geometry has the needed subsets of the probelm
-sh = dom.subset_handler()
-for e in requiredSubsets:
-    if sh.get_subset_index(e) == -1:
-        print(f"Domain does not contain subset {e}.")
-        sys.exit(1)
+
+# In[4]:
+
+
+dom = util.CreateDomain(gridName, numRefs, requiredSubsets)
+
+
+# ### Create Approximation space
+
+# In[5]:
+
 
 # Create approximation space which describes the unknowns in the equation
-fct = "u"  # name of the function
-type = "Lagrange"
-order = 1  # polynom order for lagrange
-
-approxSpace = ug4.ApproximationSpace2d(dom)
-approxSpace.add_fct(fct, type, order)
-approxSpace.init_levels()
-approxSpace.init_surfaces()  
-approxSpace.init_top_surface()
-approxSpace.print_statistic()
-
-# Create discretization for the subsets LIP and COR
-# perform the discretization on the actual elements
-# using first order finite volumes (FV1) on the 2d grid
-KDcor=1.0
-KCor=1.0
+approxSpaceDesc = dict(fct = "u", type = "Lagrange", order = 1)
+
+
+# In[6]:
+
+
+approxSpace = util.CreateApproximationSpace(dom, approxSpaceDesc)
+
+
+# ## Create a convection-diffusion-equation
+# $$\frac{\partial Ku}{\partial t} + \nabla \cdot [-DK \nabla u] = 0$$
+# Partition coefficients
+
+# In[7]:
+
+
+K={}
+K["COR"]=1e-3
+K["LIP"]=1.0
+
+
+# Diffusion coefficients
+
+# In[8]:
+
+
+D={}
+D["COR"]=1.0
+D["LIP"]=1.0
+
+
+# Creating two instances of a convection diffusion equation
+
+# In[9]:
+
 
 elemDisc={}
-# creating instance of a convection diffusion equation
-elemDisc["COR"] = cd.ConvectionDiffusionFV12d("u", "COR")
-elemDisc["COR"].set_diffusion(KDcor)
-elemDisc["COR"].set_mass_scale(KCor)
+elemDisc["COR"] = util.CreateDiffusionElemDisc("u", "COR", K["COR"], K["COR"]*D["COR"], 0.0)
+elemDisc["LIP"] = util.CreateDiffusionElemDisc("u", "LIP", K["LIP"], K["LIP"]*D["LIP"], 0.0)
 
-elemDisc["LIP"] = cd.ConvectionDiffusionFV12d("u", "LIP")
-elemDisc["LIP"].set_diffusion(KDcor)
-elemDisc["LIP"].set_mass_scale(KCor)
 
+# ### Boundary conditions
 # ug4 separates the boundary value and the discretization
 # boundary conditions can be enforced through a post-process (dirichlet).
 # To init at boundary, the value, function name from the Approximationspace
-# and the subset name are needed
+# #and the subset name are needed
+
+# In[10]:
+
+
 dirichletBND = ug4.DirichletBoundary2dCPU1()
 dirichletBND.add(1.0, "u", "TOP_SC")
 dirichletBND.add(0.0, "u", "BOTTOM_SC")
 
 
-# create the discretization object which combines all the
-# separate discretizations into one domain discretization.
+# ### Global discretization
+# Create the discretization object which combines all the separate discretizations into one domain discretization.
+
+# In[11]:
+
+
 domainDisc = ug4.DomainDiscretization2dCPU1(approxSpace)
 domainDisc.add(elemDisc["COR"])
 domainDisc.add(elemDisc["LIP"])
 domainDisc.add(dirichletBND)
 
 
-# Create Solver
-# In this case we use an LU (Lower Upper) Solver for an exact solution
+# ## Create  solver
+
+# In[12]:
+
+
 lsolver=ug4.LUCPU1()
+
+
+# In[13]:
+
+
+# import pysuperlu as slu
 # lsolver=slu.SuperLUCPU1()
 
-# Solve the transient problem
-# Use the Approximationspace to
-# create a vector of unknowns and a vector
-# which contains the right hand side
-usol = ug4.GridFunction2dCPU1(approxSpace)
 
-# Init the vector representing the unknowns with 0 to obtain
-# reproducable results
-ug4.Interpolate(0.0, usol, "u")
+# ## Solve transient problem
+# 
+# Solve the transient problem.
+# Use the Approximationspace to create a vector of unknowns and a vector which contains the right hand side
+
+# In[14]:
+
 
 # Define start time, end time and step size
 startTime = 0.0
 endTime = 1000.0
 dt = 25.0
 
-# create a time discretization with the theta-scheme
-# takes in domain and a theta
-# theta = 1.0 -> Implicit Euler, 
-# 0.0 -> Explicit Euler 
+
+# Create a time discretization (with the theta-scheme).
+# Requires domain and theta ($\theta$ = 1.0 -> implicit Euler, 
+# $\theta$ = 0.0 -> explicit Euler )
+# 
+
+# In[15]:
+
+
 timeDisc=ug4.ThetaTimeStepCPU1(domainDisc, 1.0) 
 
-# creating Time Iterator, setting the solver and step size
-timeInt = limex.LinearTimeIntegrator2dCPU1(timeDisc)
-#timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
+
+# Create time integrator (requires solver and step size).
+
+# In[16]:
+
+
+#timeInt = limex.LinearTimeIntegrator2dCPU1(timeDisc)
+timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
 timeInt.set_linear_solver(lsolver)
 timeInt.set_time_step(dt)
 
 
+# Solving the transient problem
+
+# In[17]:
+
+
+# Create grid function for solution.
+usol = ug4.GridFunction2dCPU1(approxSpace)
+
+# Init the vector representing the unknowns with 0 to obtain
+# reproducable results
+ug4.Interpolate(0.0, usol, "u")
+
 # Exporting the result to a vtu-file
 # can be visualized in paraview or with a python extension
-# We add a corresponding observer to the time integrator.
-def MyCallback(usol, step, time, dt) :
-    ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_"+str(int(step)).zfill(5)+".vtu")
-pyobserver =ug4.PythonCallbackObserver2dCPU1(MyCallback) 
-timeInt.attach_observer(pyobserver)
+ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_Initial")
+
+
+# In[18]:
+
 
-# Solve the transient problem.
 try:
     timeInt.apply(usol, endTime, usol, startTime)
 except Exception as inst:
     print(inst)
 
 
+# In[19]:
+
+
+# Write last solution.
+ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_Final")
 
-# Plot the result (if Pyvista is present)
-hasPyvista = True
+
+# ## Plotting the result using pyvista
+
+# In[20]:
+
+
+hasPyvista=True
 try:
     import pyvista
 except Exception as inst:
-    hasPyvista = False
+    print(inst)
+    hasPyvista=False
+
+
+# In[21]:
+
+
+if hasPyvista:
+    pyvista.start_xvfb()
+    result = pyvista.read('vtk/SkinDiffusion_Final.vtu')
+
+    pyvista.set_jupyter_backend('static')
+    #pyvista.set_jupyter_backend('trame')
 
-if hasPyvista :
-    result = pyvista.read('vtk/SkinDiffusion_00040.vtu')
     print("Pyvista input: ")
-    print(result)
-    result.plot(scalars="u", show_edges=True, cmap='coolwarm')
-#scalars="node_value", categories=True
+    result.plot(scalars="u", show_edges=True, cmap='jet')
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
diff --git a/content/skin/SkinDiffusionWithBinding.ipynb b/content/skin/SkinDiffusionWithBinding.ipynb
index 630eaf5..1a8aeca 100644
--- a/content/skin/SkinDiffusionWithBinding.ipynb
+++ b/content/skin/SkinDiffusionWithBinding.ipynb
@@ -4,43 +4,40 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "[<img src=\"../header.svg\">](../index.ipynd)\n",
+    "[<img src=\"../../header.svg\">](../../index.ipynd)\n",
     "\n",
     "---\n",
     "\n",
-    "# Simulation of Drug Transport across a Virtual Skin Membrane"
+    "# Drug Transport across a Virtual Skin Membrane: Binding of substances.\n",
+    "\n",
+    "In this example we solve a coupled system\n",
+    "\n",
+    "$$ \\partial_t c + \\nabla \\cdot [-D \\nabla c] + k_f \\cdot c - k_r \\cdot b = 0 $$\n",
+    "$$ \\partial_t b   - k_f \\cdot c + k_r \\cdot b = 0 $$\n",
+    "\n",
+    "The steady-state is given by $ b = \\frac{k_f}{k_r} c $\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Initialize UG4 (for 2D and standard algebra)"
+    "### Initialize UG4 "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 85,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "* Initializing: paths... done, bridge... done, plugins... done                 *\n",
-       "Loading Plugin Lua Content from /Users/anaegel/Software/ug4-git/plugins/Poroelasticity/lua/lua-include.lua\n",
-       "Loading Plugin Lua Content from /Users/anaegel/Software/ug4-git/plugins/d3f/lua/lua-include.lua\n",
-       "Loading Plugin Lua Content from /Users/anaegel/Software/ug4-git/plugins/Limex/lua/lua-include.lua\n",
-       "LIMEX-Plugin: Loading util.limex.*\n",
-       "Loading Plugin Lua Content from /Users/anaegel/Software/ug4-git/plugins/RepoTrend/lua/lua-include.lua\n",
-       "RepoTrend-Plugin: Extended package.path:./?.lua;/usr/local/share/lua/5.1/?.lua;/usr/local/share/lua/5.1/?/init.lua;/usr/local/lib/lua/5.1/?.lua;/usr/local/lib/lua/5.1/?/init.lua;/Users/anaegel/Software/ug4-git/plugins/RepoTrend/lua/?.lua\n"
-      ]
-     },
-     "execution_count": 1,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
+    "# Importing \n",
+    "import ug4py.pyugcore  as ugcore\n",
+    "import ug4py.pylimex as limex\n",
+    "import ug4py.pyconvectiondiffusion as cd\n",
+    "\n",
+    "import sys\n",
+    "sys.path.append('..')\n",
     "import modsimtools as util"
    ]
   },
@@ -53,18 +50,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 69,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 2,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
     "requiredSubsets = {\"LIP\", \"COR\", \"BOTTOM_SC\", \"TOP_SC\"}\n",
     "gridName = \"skin2d-aniso.ugx\"\n",
@@ -73,20 +61,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 70,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "Loading Domain skin2d-aniso.ugx ... done.\n",
-       "Performing integrity check on domain ... done.\n",
-       "Refining(2): 1 2 done.\n"
-      ]
-     },
-     "execution_count": 3,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Loading Domain 'skin2d-aniso.ugx'...\n",
+      "Domain loaded.\n",
+      "Refining ...\n",
+      "Refining step {0} ...\n",
+      "Refining step {1} ...\n",
+      "Refining done\n"
+     ]
     }
    ],
    "source": [
@@ -102,44 +90,41 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 71,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "Approximation space:\n",
-       "| ----------------------------------------------------------------------------------------- |\n",
-       "|  Number of DoFs (All Procs)                                                               |\n",
-       "|  Algebra: Block 1 (divide by 1 for #Index)                                                |\n",
-       "|                                                                                           |\n",
-       "|    GridLevel   |       Domain |       0: LIP |       1: COR | 2: BOTTOM_SC |    3: TOP_SC |\n",
-       "| ----------------------------------------------------------------------------------------- |\n",
-       "| (lev,    0)    |         1360 |           64 |         1216 |           40 |           40 |\n",
-       "| (lev,    1)    |         5226 |         1566 |         3504 |           78 |           78 |\n",
-       "| (lev,    2)    |        20482 |         8734 |        11440 |          154 |          154 |\n",
-       "| (lev,    0, g) |         1360 |           64 |         1216 |           40 |           40 |\n",
-       "| (lev,    1, g) |         5226 |         1566 |         3504 |           78 |           78 |\n",
-       "| (lev,    2, g) |        20482 |         8734 |        11440 |          154 |          154 |\n",
-       "| (surf, top)    |        20482 |         8734 |        11440 |          154 |          154 |\n",
-       "| ----------------------------------------------------------------------------------------- |\n",
-       "| For Legend and Options: print_statistic(\"legend\").                                        |\n",
-       "| ----------------------------------------------------------------------------------------- |\n"
-      ]
-     },
-     "execution_count": 4,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Approximation space:\n",
+      "| ----------------------------------------------------------------------------------------- |\n",
+      "|  Number of DoFs (All Procs)                                                               |\n",
+      "|  Algebra: Block 1 (divide by 1 for #Index)                                                |\n",
+      "|                                                                                           |\n",
+      "|    GridLevel   |       Domain |       0: LIP |       1: COR | 2: BOTTOM_SC |    3: TOP_SC |\n",
+      "| ----------------------------------------------------------------------------------------- |\n",
+      "| (lev,    0)    |         1360 |           64 |         1216 |           40 |           40 |\n",
+      "| (lev,    1)    |         5226 |         1566 |         3504 |           78 |           78 |\n",
+      "| (lev,    2)    |        20482 |         8734 |        11440 |          154 |          154 |\n",
+      "| (lev,    0, g) |         1360 |           64 |         1216 |           40 |           40 |\n",
+      "| (lev,    1, g) |         5226 |         1566 |         3504 |           78 |           78 |\n",
+      "| (lev,    2, g) |        20482 |         8734 |        11440 |          154 |          154 |\n",
+      "| (surf, top)    |        20482 |         8734 |        11440 |          154 |          154 |\n",
+      "| ----------------------------------------------------------------------------------------- |\n",
+      "| For Legend and Options: print_statistic(\"legend\").                                        |\n",
+      "| ----------------------------------------------------------------------------------------- |\n"
+     ]
     }
    ],
    "source": [
-    "approxSpace = ApproximationSpace(dom)\n",
-    "approxSpace:add_fct(\"c\",\"Lagrange\", 1) \n",
-    "approxSpace:add_fct(\"b\",\"Lagrange\", 1) -- NEW: Two functions are required\n",
-    "approxSpace:init_levels()\n",
-    "approxSpace:init_top_surface()\n",
+    "approxSpace = ugcore.ApproximationSpace2d(dom)\n",
+    "approxSpace.add_fct(\"c\",\"Lagrange\", 1) \n",
+    "approxSpace.add_fct(\"b\",\"Lagrange\", 1) ## NEW: Two functions are required\n",
+    "approxSpace.init_levels()\n",
+    "approxSpace.init_top_surface()\n",
     "print(\"Approximation space:\")\n",
-    "approxSpace:print_statistic()"
+    "approxSpace.print_statistic()"
    ]
   },
   {
@@ -148,34 +133,17 @@
    "source": [
     "## Create a convection-diffusion-equation\n",
     "\n",
-    "$$ \\partial_t c + \\nabla \\cdot [-D \\nabla c] + k_f \\cdot c - k_r \\cdot b = 0 $$\n",
-    "$$ \\partial_t b   - k_f \\cdot c + k_r \\cdot b = 0 $$\n",
-    "\n",
-    "with steady-state: \n",
-    "$ b = \\frac{k_f}{k_r} c $\n",
-    "\n",
     "* Define model parameters"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 72,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 5,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "D={\n",
-    "    [\"LIP\"] = 1, \n",
-    "    [\"COR\"] = 0.01, \n",
-    "}\n",
+    "# Define model parameter\n",
+    "D={ \"COR\":  0.01, \"LIP\":1.0 }\n",
     "\n",
     "kf = 2.0\n",
     "kr = 1.0"
@@ -185,36 +153,33 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "* Print characteristic times"
+    "* Print characteristic times:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 73,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "Characteristic times:\n",
-       "Binding (c->b):\t\t0.5\n",
-       "Binding (b->c):\t\t1\n",
-       "Diffusion [LIP]:\t309.76\n",
-       "Diffusion [COR]:\t30976\n"
-      ]
-     },
-     "execution_count": 6,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Characteristic times:\n",
+      "Binding (c->b):\t\t0.5\n",
+      "Binding (b->c):\t\t1.0\n",
+      "Diffusion [LIP]:\t309.76000000000005\n",
+      "Diffusion [COR]:\t30976.000000000004\n"
+     ]
     }
    ],
    "source": [
-    "L=17.6 -- characteristic length\n",
+    "L=17.6 # characteristic length\n",
     "print(\"Characteristic times:\")\n",
-    "print(\"Binding (c->b):\\t\\t\"..1.0/kf)\n",
-    "print(\"Binding (b->c):\\t\\t\"..1.0/kr)\n",
-    "print(\"Diffusion [LIP]:\\t\"..L*L/D[\"LIP\"])\n",
-    "print(\"Diffusion [COR]:\\t\"..L*L/D[\"COR\"])"
+    "print(\"Binding (c->b):\\t\\t\" + str(1.0/kf))\n",
+    "print(\"Binding (b->c):\\t\\t\" + str(1.0/kr))\n",
+    "print(\"Diffusion [LIP]:\\t\" + str(L*L/D[\"LIP\"]))\n",
+    "print(\"Diffusion [COR]:\\t\" + str(L*L/D[\"COR\"]))"
    ]
   },
   {
@@ -226,54 +191,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 74,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 7,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
+    "# Create element discretizations for free species \n",
     "elemDiscC ={}\n",
-    "\n",
-    "elemDiscC[\"LIP\"] = ConvectionDiffusion(\"c\", \"LIP\", \"fv1\")\n",
-    "elemDiscC[\"LIP\"]:set_diffusion(D[\"LIP\"])\n",
-    "elemDiscC[\"LIP\"]:set_mass_scale(1.0)\n",
-    "\n",
-    "elemDiscC[\"COR\"] = ConvectionDiffusion(\"c\", \"COR\", \"fv1\")\n",
-    "elemDiscC[\"COR\"]:set_diffusion(D[\"COR\"])\n",
-    "elemDiscC[\"COR\"]:set_mass_scale(1.0)"
+    "elemDiscC[\"COR\"] = util.CreateDiffusionElemDisc(\"c\", \"COR\", 1.0, D[\"COR\"], 0.0)\n",
+    "elemDiscC[\"LIP\"] = util.CreateDiffusionElemDisc(\"c\", \"LIP\", 1.0, D[\"LIP\"], 0.0)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 75,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 8,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
+    "# Create element discretizations for bound species\n",
     "elemDiscB ={}\n",
-    "\n",
-    "elemDiscB[\"LIP\"] = ConvectionDiffusion(\"b\", \"LIP\", \"fv1\")\n",
-    "elemDiscB[\"LIP\"]:set_diffusion(0.0)\n",
-    "elemDiscB[\"LIP\"]:set_mass_scale(1.0)\n",
-    "\n",
-    "elemDiscB[\"COR\"] = ConvectionDiffusion(\"b\", \"COR\", \"fv1\")\n",
-    "elemDiscB[\"COR\"]:set_diffusion(0.0)\n",
-    "elemDiscB[\"COR\"]:set_mass_scale(1.0)"
+    "elemDiscB[\"COR\"] = util.CreateDiffusionElemDisc(\"b\", \"COR\", 1.0, 0.0, 0.0)\n",
+    "elemDiscB[\"LIP\"] = util.CreateDiffusionElemDisc(\"b\", \"LIP\", 1.0, 0.0, 0.0)"
    ]
   },
   {
@@ -285,24 +222,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 76,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 9,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "myReactionC = kf*elemDiscC[\"COR\"]:value()-kr*elemDiscB[\"COR\"]:value()\n",
-    "myReactionB = -kf*elemDiscC[\"COR\"]:value()+kr*elemDiscB[\"COR\"]:value()\n",
+    "#myReactionC = kf*elemDiscC[\"COR\"]:value()-kr*elemDiscB[\"COR\"]:value()\n",
+    "#myReactionB = -kf*elemDiscC[\"COR\"]:value()+kr*elemDiscB[\"COR\"]:value()\n",
+    "\n",
+    "def MyFunction (b,c): return kf*c - kr*b\n",
+    "def MyFunction_b (b,c): return kf\n",
+    "def MyFunction_c (b,c): return -kr\n",
+    "\n",
+    "pyReactionC=ugcore.PythonUserFunction2d(MyFunction, 2)\n",
+    "pyReactionC.set_input_and_deriv(0,elemDiscC[\"COR\"].value(), MyFunction_b)\n",
+    "pyReactionC.set_input_and_deriv(1,elemDiscB[\"COR\"].value(), MyFunction_c)\n",
     "\n",
-    "elemDiscC[\"COR\"]:set_reaction(myReactionC)\n",
-    "elemDiscB[\"COR\"]:set_reaction(myReactionB)"
+    "pyReactionB =  ugcore.ScaleAddLinkerNumber2d()\n",
+    "pyReactionB.add(-1.0, pyReactionC)\n",
+    "\n",
+    "elemDiscC[\"COR\"].set_reaction(pyReactionC)\n",
+    "elemDiscB[\"COR\"].set_reaction(pyReactionB)"
    ]
   },
   {
@@ -314,24 +253,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 77,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 10,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "dirichletBnd = DirichletBoundary()\n",
-    "dirichletBnd:add(1.0, \"c\", \"TOP_SC\")\n",
-    "dirichletBnd:add(0.0, \"c\", \"BOTTOM_SC\")\n",
-    "dirichletBnd:add(0.0, \"b\", \"TOP_SC\") -- NEW\n",
-    "dirichletBnd:add(0.0, \"b\", \"BOTTOM_SC\")"
+    "dirichletBnd = ugcore.DirichletBoundary2dCPU1()\n",
+    "dirichletBnd.add(1.0, \"c\", \"TOP_SC\")\n",
+    "dirichletBnd.add(0.0, \"c\", \"BOTTOM_SC\")\n",
+    "dirichletBnd.add(0.0, \"b\", \"TOP_SC\") #NEW\n",
+    "dirichletBnd.add(0.0, \"b\", \"BOTTOM_SC\")"
    ]
   },
   {
@@ -343,52 +273,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 11,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "domainDisc = DomainDiscretization(approxSpace)\n",
-    "domainDisc:add(elemDiscC[\"LIP\"]) \n",
-    "domainDisc:add(elemDiscC[\"COR\"])\n",
-    "\n",
-    "-- NEW: More discs\n",
-    "domainDisc:add(elemDiscB[\"LIP\"])\n",
-    "domainDisc:add(elemDiscB[\"COR\"])\n",
-    "domainDisc:add(dirichletBnd)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Create  solver"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 78,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 12,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "lsolver=SuperLU()"
+    "try:\n",
+    "    domainDisc = ugcore.DomainDiscretization2dCPU1(approxSpace)\n",
+    "    domainDisc.add(dirichletBnd)\n",
+    "    domainDisc.add(elemDiscC[\"LIP\"]) \n",
+    "    domainDisc.add(elemDiscC[\"COR\"])\n",
+    "    domainDisc.add(elemDiscB[\"LIP\"])\n",
+    "    domainDisc.add(elemDiscB[\"COR\"])\n",
+    "except RuntimeError as inst:    \n",
+    "    print(type(inst))    # the exception type\n",
+    "    print(inst.args)     # arguments stored in .args\n",
+    "    print(inst)  "
    ]
   },
   {
@@ -401,18 +300,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 79,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 13,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
     "startTime = 0.0\n",
     "endTime = 500.0\n",
@@ -428,554 +318,109 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 80,
    "metadata": {
     "scrolled": true
    },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 14,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "u = GridFunction(approxSpace)\n",
-    "u:set(0.0)"
+    "usol = ugcore.GridFunction2dCPU1(approxSpace)\n",
+    "ugcore.Interpolate(0.0, usol, \"c\")\n",
+    "ugcore.Interpolate(0.0, usol, \"b\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "* Solve problem (using lua util):"
+    "* Create solver"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 81,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "SolveLinearTimeProblem, Linear Solver setup:\n",
-       "SuperLU\n",
-       ">> Writing start values\n",
-       "++++++ TIMESTEP 1 BEGIN (current time: 0) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ Assembling Matrix/Rhs for step size 5\n",
-       "++++++ TIMESTEP 1 END   (current time: 5) ++++++\n",
-       "++++++ TIMESTEP 2 BEGIN (current time: 5) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 2 END   (current time: 10) ++++++\n",
-       "++++++ TIMESTEP 3 BEGIN (current time: 10) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 3 END   (current time: 15) ++++++\n",
-       "++++++ TIMESTEP 4 BEGIN (current time: 15) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 4 END   (current time: 20) ++++++\n",
-       "++++++ TIMESTEP 5 BEGIN (current time: 20) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 5 END   (current time: 25) ++++++\n",
-       "++++++ TIMESTEP 6 BEGIN (current time: 25) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 6 END   (current time: 30) ++++++\n",
-       "++++++ TIMESTEP 7 BEGIN (current time: 30) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 7 END   (current time: 35) ++++++\n",
-       "++++++ TIMESTEP 8 BEGIN (current time: 35) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 8 END   (current time: 40) ++++++\n",
-       "++++++ TIMESTEP 9 BEGIN (current time: 40) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 9 END   (current time: 45) ++++++\n",
-       "++++++ TIMESTEP 10 BEGIN (current time: 45) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 10 END   (current time: 50) ++++++\n",
-       "++++++ TIMESTEP 11 BEGIN (current time: 50) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 11 END   (current time: 55) ++++++\n",
-       "++++++ TIMESTEP 12 BEGIN (current time: 55) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 12 END   (current time: 60) ++++++\n",
-       "++++++ TIMESTEP 13 BEGIN (current time: 60) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 13 END   (current time: 65) ++++++\n",
-       "++++++ TIMESTEP 14 BEGIN (current time: 65) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 14 END   (current time: 70) ++++++\n",
-       "++++++ TIMESTEP 15 BEGIN (current time: 70) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 15 END   (current time: 75) ++++++\n",
-       "++++++ TIMESTEP 16 BEGIN (current time: 75) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 16 END   (current time: 80) ++++++\n",
-       "++++++ TIMESTEP 17 BEGIN (current time: 80) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 17 END   (current time: 85) ++++++\n",
-       "++++++ TIMESTEP 18 BEGIN (current time: 85) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 18 END   (current time: 90) ++++++\n",
-       "++++++ TIMESTEP 19 BEGIN (current time: 90) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 19 END   (current time: 95) ++++++\n",
-       "++++++ TIMESTEP 20 BEGIN (current time: 95) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 20 END   (current time: 100) ++++++\n",
-       "++++++ TIMESTEP 21 BEGIN (current time: 100) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 21 END   (current time: 105) ++++++\n",
-       "++++++ TIMESTEP 22 BEGIN (current time: 105) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 22 END   (current time: 110) ++++++\n",
-       "++++++ TIMESTEP 23 BEGIN (current time: 110) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 23 END   (current time: 115) ++++++\n",
-       "++++++ TIMESTEP 24 BEGIN (current time: 115) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 24 END   (current time: 120) ++++++\n",
-       "++++++ TIMESTEP 25 BEGIN (current time: 120) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 25 END   (current time: 125) ++++++\n",
-       "++++++ TIMESTEP 26 BEGIN (current time: 125) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 26 END   (current time: 130) ++++++\n",
-       "++++++ TIMESTEP 27 BEGIN (current time: 130) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 27 END   (current time: 135) ++++++\n",
-       "++++++ TIMESTEP 28 BEGIN (current time: 135) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 28 END   (current time: 140) ++++++\n",
-       "++++++ TIMESTEP 29 BEGIN (current time: 140) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 29 END   (current time: 145) ++++++\n",
-       "++++++ TIMESTEP 30 BEGIN (current time: 145) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 30 END   (current time: 150) ++++++\n",
-       "++++++ TIMESTEP 31 BEGIN (current time: 150) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 31 END   (current time: 155) ++++++\n",
-       "++++++ TIMESTEP 32 BEGIN (current time: 155) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 32 END   (current time: 160) ++++++\n",
-       "++++++ TIMESTEP 33 BEGIN (current time: 160) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 33 END   (current time: 165) ++++++\n",
-       "++++++ TIMESTEP 34 BEGIN (current time: 165) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 34 END   (current time: 170) ++++++\n",
-       "++++++ TIMESTEP 35 BEGIN (current time: 170) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 35 END   (current time: 175) ++++++\n",
-       "++++++ TIMESTEP 36 BEGIN (current time: 175) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 36 END   (current time: 180) ++++++\n",
-       "++++++ TIMESTEP 37 BEGIN (current time: 180) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 37 END   (current time: 185) ++++++\n",
-       "++++++ TIMESTEP 38 BEGIN (current time: 185) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 38 END   (current time: 190) ++++++\n",
-       "++++++ TIMESTEP 39 BEGIN (current time: 190) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 39 END   (current time: 195) ++++++\n",
-       "++++++ TIMESTEP 40 BEGIN (current time: 195) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 40 END   (current time: 200) ++++++\n",
-       "++++++ TIMESTEP 41 BEGIN (current time: 200) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 41 END   (current time: 205) ++++++\n",
-       "++++++ TIMESTEP 42 BEGIN (current time: 205) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 42 END   (current time: 210) ++++++\n",
-       "++++++ TIMESTEP 43 BEGIN (current time: 210) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 43 END   (current time: 215) ++++++\n",
-       "++++++ TIMESTEP 44 BEGIN (current time: 215) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 44 END   (current time: 220) ++++++\n",
-       "++++++ TIMESTEP 45 BEGIN (current time: 220) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 45 END   (current time: 225) ++++++\n",
-       "++++++ TIMESTEP 46 BEGIN (current time: 225) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 46 END   (current time: 230) ++++++\n",
-       "++++++ TIMESTEP 47 BEGIN (current time: 230) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 47 END   (current time: 235) ++++++\n",
-       "++++++ TIMESTEP 48 BEGIN (current time: 235) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 48 END   (current time: 240) ++++++\n",
-       "++++++ TIMESTEP 49 BEGIN (current time: 240) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 49 END   (current time: 245) ++++++\n",
-       "++++++ TIMESTEP 50 BEGIN (current time: 245) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 50 END   (current time: 250) ++++++\n",
-       "++++++ TIMESTEP 51 BEGIN (current time: 250) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 51 END   (current time: 255) ++++++\n",
-       "++++++ TIMESTEP 52 BEGIN (current time: 255) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 52 END   (current time: 260) ++++++\n",
-       "++++++ TIMESTEP 53 BEGIN (current time: 260) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 53 END   (current time: 265) ++++++\n",
-       "++++++ TIMESTEP 54 BEGIN (current time: 265) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 54 END   (current time: 270) ++++++\n",
-       "++++++ TIMESTEP 55 BEGIN (current time: 270) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 55 END   (current time: 275) ++++++\n",
-       "++++++ TIMESTEP 56 BEGIN (current time: 275) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 56 END   (current time: 280) ++++++\n",
-       "++++++ TIMESTEP 57 BEGIN (current time: 280) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 57 END   (current time: 285) ++++++\n",
-       "++++++ TIMESTEP 58 BEGIN (current time: 285) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 58 END   (current time: 290) ++++++\n",
-       "++++++ TIMESTEP 59 BEGIN (current time: 290) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 59 END   (current time: 295) ++++++\n",
-       "++++++ TIMESTEP 60 BEGIN (current time: 295) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 60 END   (current time: 300) ++++++\n",
-       "++++++ TIMESTEP 61 BEGIN (current time: 300) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 61 END   (current time: 305) ++++++\n",
-       "++++++ TIMESTEP 62 BEGIN (current time: 305) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 62 END   (current time: 310) ++++++\n",
-       "++++++ TIMESTEP 63 BEGIN (current time: 310) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 63 END   (current time: 315) ++++++\n",
-       "++++++ TIMESTEP 64 BEGIN (current time: 315) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 64 END   (current time: 320) ++++++\n",
-       "++++++ TIMESTEP 65 BEGIN (current time: 320) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 65 END   (current time: 325) ++++++\n",
-       "++++++ TIMESTEP 66 BEGIN (current time: 325) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 66 END   (current time: 330) ++++++\n",
-       "++++++ TIMESTEP 67 BEGIN (current time: 330) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 67 END   (current time: 335) ++++++\n",
-       "++++++ TIMESTEP 68 BEGIN (current time: 335) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 68 END   (current time: 340) ++++++\n",
-       "++++++ TIMESTEP 69 BEGIN (current time: 340) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 69 END   (current time: 345) ++++++\n",
-       "++++++ TIMESTEP 70 BEGIN (current time: 345) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 70 END   (current time: 350) ++++++\n",
-       "++++++ TIMESTEP 71 BEGIN (current time: 350) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 71 END   (current time: 355) ++++++\n",
-       "++++++ TIMESTEP 72 BEGIN (current time: 355) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 72 END   (current time: 360) ++++++\n",
-       "++++++ TIMESTEP 73 BEGIN (current time: 360) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 73 END   (current time: 365) ++++++\n",
-       "++++++ TIMESTEP 74 BEGIN (current time: 365) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 74 END   (current time: 370) ++++++\n",
-       "++++++ TIMESTEP 75 BEGIN (current time: 370) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 75 END   (current time: 375) ++++++\n",
-       "++++++ TIMESTEP 76 BEGIN (current time: 375) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 76 END   (current time: 380) ++++++\n",
-       "++++++ TIMESTEP 77 BEGIN (current time: 380) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 77 END   (current time: 385) ++++++\n",
-       "++++++ TIMESTEP 78 BEGIN (current time: 385) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 78 END   (current time: 390) ++++++\n",
-       "++++++ TIMESTEP 79 BEGIN (current time: 390) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 79 END   (current time: 395) ++++++\n",
-       "++++++ TIMESTEP 80 BEGIN (current time: 395) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 80 END   (current time: 400) ++++++\n",
-       "++++++ TIMESTEP 81 BEGIN (current time: 400) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 81 END   (current time: 405) ++++++\n",
-       "++++++ TIMESTEP 82 BEGIN (current time: 405) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 82 END   (current time: 410) ++++++\n",
-       "++++++ TIMESTEP 83 BEGIN (current time: 410) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 83 END   (current time: 415) ++++++\n",
-       "++++++ TIMESTEP 84 BEGIN (current time: 415) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 84 END   (current time: 420) ++++++\n",
-       "++++++ TIMESTEP 85 BEGIN (current time: 420) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 85 END   (current time: 425) ++++++\n",
-       "++++++ TIMESTEP 86 BEGIN (current time: 425) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 86 END   (current time: 430) ++++++\n",
-       "++++++ TIMESTEP 87 BEGIN (current time: 430) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 87 END   (current time: 435) ++++++\n",
-       "++++++ TIMESTEP 88 BEGIN (current time: 435) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 88 END   (current time: 440) ++++++\n",
-       "++++++ TIMESTEP 89 BEGIN (current time: 440) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 89 END   (current time: 445) ++++++\n",
-       "++++++ TIMESTEP 90 BEGIN (current time: 445) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 90 END   (current time: 450) ++++++\n",
-       "++++++ TIMESTEP 91 BEGIN (current time: 450) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 91 END   (current time: 455) ++++++\n",
-       "++++++ TIMESTEP 92 BEGIN (current time: 455) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 92 END   (current time: 460) ++++++\n",
-       "++++++ TIMESTEP 93 BEGIN (current time: 460) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 93 END   (current time: 465) ++++++\n",
-       "++++++ TIMESTEP 94 BEGIN (current time: 465) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 94 END   (current time: 470) ++++++\n",
-       "++++++ TIMESTEP 95 BEGIN (current time: 470) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 95 END   (current time: 475) ++++++\n",
-       "++++++ TIMESTEP 96 BEGIN (current time: 475) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 96 END   (current time: 480) ++++++\n",
-       "++++++ TIMESTEP 97 BEGIN (current time: 480) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 97 END   (current time: 485) ++++++\n",
-       "++++++ TIMESTEP 98 BEGIN (current time: 485) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 98 END   (current time: 490) ++++++\n",
-       "++++++ TIMESTEP 99 BEGIN (current time: 490) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 99 END   (current time: 495) ++++++\n",
-       "++++++ TIMESTEP 100 BEGIN (current time: 495) ++++++\n",
-       "++++++ Time step size: 5\n",
-       "++++++ TIMESTEP 100 END   (current time: 500) ++++++\n"
-      ]
-     },
-     "execution_count": 15,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "util.SolveLinearTimeProblem(u, domainDisc, lsolver, VTKOutput(), \"vtk/bindingLo\",\n",
-    "                            \"ImplEuler\", 1, startTime, endTime, dt); "
+    "lsolver=ugcore.LUCPU1()"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Alternative\n",
     "* Use C++-Objects"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 82,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 16,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "local timeDisc=ThetaTimeStep(domainDisc, 1.0)\n",
-    "timeInt = ConstStepLinearTimeIntegrator(timeDisc)\n",
-    "timeInt:set_linear_solver(lsolver)\n",
-    "timeInt:set_time_step(dt)"
+    "timeDisc=ugcore.ThetaTimeStepCPU1(domainDisc, 1.0)\n",
+    "timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)\n",
+    "timeInt.set_linear_solver(lsolver)\n",
+    "timeInt.set_time_step(dt)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 83,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": []
-     },
-     "execution_count": 17,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "local vtkObserver = VTKOutputObserver(\"vtk/BindingSol.vtk\", VTKOutput())\n",
-    "timeInt:attach_finalize_observer(vtkObserver)"
+    "vtkObserver = limex.VTKOutputObserver2dCPU1(\"vtk/BindingSol.vtk\", ugcore.VTKOutput2d())\n",
+    "timeInt.attach_finalize_observer(vtkObserver)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 84,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "+++ Integrating: [\t0\t, \t500\t] with dt=\t5(100 iters)\n",
-       "+++ Const timestep +++1(t=0, dt=5)\n",
-       "+++ Assemble (t=0, dt=5)\n",
-       "+++ Const timestep +++2(t=5, dt=5)\n",
-       "+++ Const timestep +++3(t=10, dt=5)\n",
-       "+++ Const timestep +++4(t=15, dt=5)\n",
-       "+++ Const timestep +++5(t=20, dt=5)\n",
-       "+++ Const timestep +++6(t=25, dt=5)\n",
-       "+++ Const timestep +++7(t=30, dt=5)\n",
-       "+++ Const timestep +++8(t=35, dt=5)\n",
-       "+++ Const timestep +++9(t=40, dt=5)\n",
-       "+++ Const timestep +++10(t=45, dt=5)\n",
-       "+++ Const timestep +++11(t=50, dt=5)\n",
-       "+++ Const timestep +++12(t=55, dt=5)\n",
-       "+++ Const timestep +++13(t=60, dt=5)\n",
-       "+++ Const timestep +++14(t=65, dt=5)\n",
-       "+++ Const timestep +++15(t=70, dt=5)\n",
-       "+++ Const timestep +++16(t=75, dt=5)\n",
-       "+++ Const timestep +++17(t=80, dt=5)\n",
-       "+++ Const timestep +++18(t=85, dt=5)\n",
-       "+++ Const timestep +++19(t=90, dt=5)\n",
-       "+++ Const timestep +++20(t=95, dt=5)\n",
-       "+++ Const timestep +++21(t=100, dt=5)\n",
-       "+++ Const timestep +++22(t=105, dt=5)\n",
-       "+++ Const timestep +++23(t=110, dt=5)\n",
-       "+++ Const timestep +++24(t=115, dt=5)\n",
-       "+++ Const timestep +++25(t=120, dt=5)\n",
-       "+++ Const timestep +++26(t=125, dt=5)\n",
-       "+++ Const timestep +++27(t=130, dt=5)\n",
-       "+++ Const timestep +++28(t=135, dt=5)\n",
-       "+++ Const timestep +++29(t=140, dt=5)\n",
-       "+++ Const timestep +++30(t=145, dt=5)\n",
-       "+++ Const timestep +++31(t=150, dt=5)\n",
-       "+++ Const timestep +++32(t=155, dt=5)\n",
-       "+++ Const timestep +++33(t=160, dt=5)\n",
-       "+++ Const timestep +++34(t=165, dt=5)\n",
-       "+++ Const timestep +++35(t=170, dt=5)\n",
-       "+++ Const timestep +++36(t=175, dt=5)\n",
-       "+++ Const timestep +++37(t=180, dt=5)\n",
-       "+++ Const timestep +++38(t=185, dt=5)\n",
-       "+++ Const timestep +++39(t=190, dt=5)\n",
-       "+++ Const timestep +++40(t=195, dt=5)\n",
-       "+++ Const timestep +++41(t=200, dt=5)\n",
-       "+++ Const timestep +++42(t=205, dt=5)\n",
-       "+++ Const timestep +++43(t=210, dt=5)\n",
-       "+++ Const timestep +++44(t=215, dt=5)\n",
-       "+++ Const timestep +++45(t=220, dt=5)\n",
-       "+++ Const timestep +++46(t=225, dt=5)\n",
-       "+++ Const timestep +++47(t=230, dt=5)\n",
-       "+++ Const timestep +++48(t=235, dt=5)\n",
-       "+++ Const timestep +++49(t=240, dt=5)\n",
-       "+++ Const timestep +++50(t=245, dt=5)\n",
-       "+++ Const timestep +++51(t=250, dt=5)\n",
-       "+++ Const timestep +++52(t=255, dt=5)\n",
-       "+++ Const timestep +++53(t=260, dt=5)\n",
-       "+++ Const timestep +++54(t=265, dt=5)\n",
-       "+++ Const timestep +++55(t=270, dt=5)\n",
-       "+++ Const timestep +++56(t=275, dt=5)\n",
-       "+++ Const timestep +++57(t=280, dt=5)\n",
-       "+++ Const timestep +++58(t=285, dt=5)\n",
-       "+++ Const timestep +++59(t=290, dt=5)\n",
-       "+++ Const timestep +++60(t=295, dt=5)\n",
-       "+++ Const timestep +++61(t=300, dt=5)\n",
-       "+++ Const timestep +++62(t=305, dt=5)\n",
-       "+++ Const timestep +++63(t=310, dt=5)\n",
-       "+++ Const timestep +++64(t=315, dt=5)\n",
-       "+++ Const timestep +++65(t=320, dt=5)\n",
-       "+++ Const timestep +++66(t=325, dt=5)\n",
-       "+++ Const timestep +++67(t=330, dt=5)\n",
-       "+++ Const timestep +++68(t=335, dt=5)\n",
-       "+++ Const timestep +++69(t=340, dt=5)\n",
-       "+++ Const timestep +++70(t=345, dt=5)\n",
-       "+++ Const timestep +++71(t=350, dt=5)\n",
-       "+++ Const timestep +++72(t=355, dt=5)\n",
-       "+++ Const timestep +++73(t=360, dt=5)\n",
-       "+++ Const timestep +++74(t=365, dt=5)\n",
-       "+++ Const timestep +++75(t=370, dt=5)\n",
-       "+++ Const timestep +++76(t=375, dt=5)\n",
-       "+++ Const timestep +++77(t=380, dt=5)\n",
-       "+++ Const timestep +++78(t=385, dt=5)\n",
-       "+++ Const timestep +++79(t=390, dt=5)\n",
-       "+++ Const timestep +++80(t=395, dt=5)\n",
-       "+++ Const timestep +++81(t=400, dt=5)\n",
-       "+++ Const timestep +++82(t=405, dt=5)\n",
-       "+++ Const timestep +++83(t=410, dt=5)\n",
-       "+++ Const timestep +++84(t=415, dt=5)\n",
-       "+++ Const timestep +++85(t=420, dt=5)\n",
-       "+++ Const timestep +++86(t=425, dt=5)\n",
-       "+++ Const timestep +++87(t=430, dt=5)\n",
-       "+++ Const timestep +++88(t=435, dt=5)\n",
-       "+++ Const timestep +++89(t=440, dt=5)\n",
-       "+++ Const timestep +++90(t=445, dt=5)\n",
-       "+++ Const timestep +++91(t=450, dt=5)\n",
-       "+++ Const timestep +++92(t=455, dt=5)\n",
-       "+++ Const timestep +++93(t=460, dt=5)\n",
-       "+++ Const timestep +++94(t=465, dt=5)\n",
-       "+++ Const timestep +++95(t=470, dt=5)\n",
-       "+++ Const timestep +++96(t=475, dt=5)\n",
-       "+++ Const timestep +++97(t=480, dt=5)\n",
-       "+++ Const timestep +++98(t=485, dt=5)\n",
-       "+++ Const timestep +++99(t=490, dt=5)\n",
-       "+++ Const timestep +++100(t=495, dt=5)\n"
-      ]
-     },
-     "execution_count": 18,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+++ Integrating: [\t0\t, \t500\t] with dt=\t5(100 iters)\n",
+      "+++ Const timestep +++1(t=0, dt=5)\n",
+      "<class 'UnicodeDecodeError'>\n",
+      "('utf-8', b'\\x01\\x00\\x00\\x00\\x00\\x00\\x00\\x00\\xd8U\\xf8\\x02\\x01\\x00\\x00\\x00\\xfe\\x04\\x00\\x00\\x00\\x00\\x00\\x00\\xff\\xff\\xff\\xff\\xff\\xff\\xff\\xffds/py-ba\\x01\\x00\\x00\\x00\\x00\\x00\\x00\\x00ls/ugcore/ugbase/lib_disc/spatial_disc/user_data/data_evaluator.cpp:191 : Unknown std::exception: TypeError: Unregistered type : std::__1::vector<double, std::__1::allocator<double>>\\n/Users/runner/work/py-basic-wheels/py-basic-wheels/ugcore/ugbase/lib_disc/spatial_disc/user_data/data_evaluator.cpp:191 : DataEvaluatorBase::prep_elem: Cannot compute data for Export or Linker.\\n/Users/runner/work/py-basic-wheels/py-basic-wheels/ugcore/cmake/../../ugcore/ugbase/lib_disc/spatial_disc/elem_disc/elem_disc_assemble_util.h:1117 : (instationary) AssembleLinear: Cannot prepare element.\\n/Users/runner/work/py-basic-wheels/py-basic-wheels/ugcore/cmake/../../ugcore/ugbase/lib_disc/spatial_disc/elem_disc/elem_disc_assemble_util.h:1225 : (instationary) AssembleLinear: Cannot create Data Evaluator.\\n/Users/runner/work/py-basic-wheels/py-basic-wheels/ugcore/cmake/../../ugcore/ugbase/lib_disc/spatial_disc/domain_disc_impl.h:1831 : DomainDiscretization::assemble_linear (instationary): Assembling of elements of Dimension 2 in  subset 1 failed.\\n/Users/runner/work/py-basic-wheels/py-basic-wheels/ugcore/ugbase/bridge/disc_bridges/../../lib_disc/time_disc/theta_time_step_impl.h:200 : MultiStepTimeDiscretization: Cannot assemble jacobian.\\n', 8, 9, 'invalid continuation byte')\n",
+      "'utf-8' codec can't decode byte 0xd8 in position 8: invalid continuation byte\n",
+      "+++ Assemble (t=0, dt=5)\n"
+     ]
     }
    ],
    "source": [
-    "timeInt:apply(u, endTime, u, startTime)"
+    "try:\n",
+    "    timeInt.apply(usol, endTime, usol, startTime)\n",
+    "except Exception as inst:    \n",
+    "    print(type(inst))    # the exception type\n",
+    "    print(inst.args)     # arguments stored in .args\n",
+    "    print(inst)  "
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "UG4-Lua 0.1.0",
-   "language": "lua",
-   "name": "ug4_kernel"
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
   },
   "language_info": {
-   "codemirror_mode": "text/x-lua",
-   "file_extension": ".lua",
-   "mimetype": "text/x-lua",
-   "name": "lua",
-   "version": "1.0"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/content/skin/SkinDiffusionWithBinding.py b/content/skin/SkinDiffusionWithBinding.py
new file mode 100644
index 0000000..014dbe2
--- /dev/null
+++ b/content/skin/SkinDiffusionWithBinding.py
@@ -0,0 +1,219 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# [<img src="../../header.svg">](../../index.ipynd)
+# 
+# ---
+# 
+# # Drug Transport across a Virtual Skin Membrane: Binding of substances.
+# 
+# In this example we solve a coupled system
+# 
+# $$ \partial_t c + \nabla \cdot [-D \nabla c] + k_f \cdot c - k_r \cdot b = 0 $$
+# $$ \partial_t b   - k_f \cdot c + k_r \cdot b = 0 $$
+# 
+# The steady-state is given by $ b = \frac{k_f}{k_r} c $
+# 
+
+# ### Initialize UG4 
+
+# In[85]:
+
+
+# Importing 
+import ug4py.pyugcore  as ugcore
+import ug4py.pylimex as limex
+import ug4py.pyconvectiondiffusion as cd
+
+import sys
+sys.path.append('..')
+import modsimtools as util
+
+
+# ### Create Domain
+
+# In[69]:
+
+
+requiredSubsets = {"LIP", "COR", "BOTTOM_SC", "TOP_SC"}
+gridName = "skin2d-aniso.ugx"
+numRefs = 2
+
+
+# In[70]:
+
+
+dom = util.CreateDomain(gridName, numRefs, requiredSubsets)
+
+
+# ### Create Approximation space
+
+# In[71]:
+
+
+approxSpace = ugcore.ApproximationSpace2d(dom)
+approxSpace.add_fct("c","Lagrange", 1) 
+approxSpace.add_fct("b","Lagrange", 1) ## NEW: Two functions are required
+approxSpace.init_levels()
+approxSpace.init_top_surface()
+print("Approximation space:")
+approxSpace.print_statistic()
+
+
+# ## Create a convection-diffusion-equation
+# 
+# * Define model parameters
+
+# In[72]:
+
+
+# Define model parameter
+D={ "COR":  0.01, "LIP":1.0 }
+
+kf = 2.0
+kr = 1.0
+
+
+# * Print characteristic times:
+
+# In[73]:
+
+
+L=17.6 # characteristic length
+print("Characteristic times:")
+print("Binding (c->b):\t\t" + str(1.0/kf))
+print("Binding (b->c):\t\t" + str(1.0/kr))
+print("Diffusion [LIP]:\t" + str(L*L/D["LIP"]))
+print("Diffusion [COR]:\t" + str(L*L/D["COR"]))
+
+
+# * The system requires **separate** element discs **for each species**:
+
+# In[74]:
+
+
+# Create element discretizations for free species 
+elemDiscC ={}
+elemDiscC["COR"] = util.CreateDiffusionElemDisc("c", "COR", 1.0, D["COR"], 0.0)
+elemDiscC["LIP"] = util.CreateDiffusionElemDisc("c", "LIP", 1.0, D["LIP"], 0.0)
+
+
+# In[75]:
+
+
+# Create element discretizations for bound species
+elemDiscB ={}
+elemDiscB["COR"] = util.CreateDiffusionElemDisc("b", "COR", 1.0, 0.0, 0.0)
+elemDiscB["LIP"] = util.CreateDiffusionElemDisc("b", "LIP", 1.0, 0.0, 0.0)
+
+
+# * Add corneocyte binding, i.e., reactions in corneocytes
+
+# In[76]:
+
+
+#myReactionC = kf*elemDiscC["COR"]:value()-kr*elemDiscB["COR"]:value()
+#myReactionB = -kf*elemDiscC["COR"]:value()+kr*elemDiscB["COR"]:value()
+
+def MyFunction (b,c): return kf*c - kr*b
+def MyFunction_b (b,c): return kf
+def MyFunction_c (b,c): return -kr
+
+pyReactionC=ugcore.PythonUserFunction2d(MyFunction, 2)
+pyReactionC.set_input_and_deriv(0,elemDiscC["COR"].value(), MyFunction_b)
+pyReactionC.set_input_and_deriv(1,elemDiscB["COR"].value(), MyFunction_c)
+
+pyReactionB =  ugcore.ScaleAddLinkerNumber2d()
+pyReactionB.add(-1.0, pyReactionC)
+
+elemDiscC["COR"].set_reaction(pyReactionC)
+elemDiscB["COR"].set_reaction(pyReactionB)
+
+
+# * Boundary conditions:
+
+# In[77]:
+
+
+dirichletBnd = ugcore.DirichletBoundary2dCPU1()
+dirichletBnd.add(1.0, "c", "TOP_SC")
+dirichletBnd.add(0.0, "c", "BOTTOM_SC")
+dirichletBnd.add(0.0, "b", "TOP_SC") #NEW
+dirichletBnd.add(0.0, "b", "BOTTOM_SC")
+
+
+# * Summarize everything in domain discretization:
+
+# In[78]:
+
+
+try:
+    domainDisc = ugcore.DomainDiscretization2dCPU1(approxSpace)
+    domainDisc.add(dirichletBnd)
+    domainDisc.add(elemDiscC["LIP"]) 
+    domainDisc.add(elemDiscC["COR"])
+    domainDisc.add(elemDiscB["LIP"])
+    domainDisc.add(elemDiscB["COR"])
+except RuntimeError as inst:    
+    print(type(inst))    # the exception type
+    print(inst.args)     # arguments stored in .args
+    print(inst)  
+
+
+# ## Solve transient problem
+# * Define parameters for time stepping:
+
+# In[79]:
+
+
+startTime = 0.0
+endTime = 500.0
+dt=endTime/100.0
+
+
+# * Define grid function:
+
+# In[80]:
+
+
+usol = ugcore.GridFunction2dCPU1(approxSpace)
+ugcore.Interpolate(0.0, usol, "c")
+ugcore.Interpolate(0.0, usol, "b")
+
+
+# * Create solver
+
+# In[81]:
+
+
+lsolver=ugcore.LUCPU1()
+
+
+# * Use C++-Objects
+
+# In[82]:
+
+
+timeDisc=ugcore.ThetaTimeStepCPU1(domainDisc, 1.0)
+timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
+timeInt.set_linear_solver(lsolver)
+timeInt.set_time_step(dt)
+
+
+# In[83]:
+
+
+vtkObserver = limex.VTKOutputObserver2dCPU1("vtk/BindingSol.vtk", ugcore.VTKOutput2d())
+timeInt.attach_finalize_observer(vtkObserver)
+
+
+# In[84]:
+
+
+try:
+    timeInt.apply(usol, endTime, usol, startTime)
+except Exception as inst:    
+    print(type(inst))    # the exception type
+    print(inst.args)     # arguments stored in .args
+    print(inst)  
+
diff --git a/content/skin/SkinDiffusionWithLagtime.ipynb b/content/skin/SkinDiffusionWithLagtime.ipynb
index bf851b8..40eaae7 100644
--- a/content/skin/SkinDiffusionWithLagtime.ipynb
+++ b/content/skin/SkinDiffusionWithLagtime.ipynb
@@ -9,7 +9,7 @@
     "---\n",
     "\n",
     "# Drug Transport across a Virtual Skin Membrane:  Extended version\n",
-    "[Previous Version](SkinDiffusion.ipybnd)"
+    "We know the basic steps from the [previous version](SkinDiffusion.ipybnd)."
    ]
   },
   {
@@ -23,16 +23,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 27,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import modsimtools as util\n",
-    "\n",
-    "import ug4py as ug4\n",
-    "import pylimex as limex\n",
-    "import pyconvectiondiffusion as cd\n",
-    "import pysuperlu as slu"
+    "import ug4py.pyugcore  as ugcore\n",
+    "import ug4py.pylimex as limex\n",
+    "import ug4py.pyconvectiondiffusion as cd\n",
+    "# import pysuperlu as slu\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Importing\n",
+    "import sys\n",
+    "sys.path.append('..')\n",
+    "import modsimtools as util"
    ]
   },
   {
@@ -44,7 +54,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 29,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -55,18 +65,18 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Loading Domain {gridName}...\n",
+      "Loading Domain 'skin2d-aniso.ugx'...\n",
       "Domain loaded.\n",
       "Refining ...\n",
-      "Refining step {i} ...\n",
-      "Refining step {i} ...\n",
+      "Refining step {0} ...\n",
+      "Refining step {1} ...\n",
       "Refining done\n"
      ]
     }
@@ -84,7 +94,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [
     {
@@ -125,7 +135,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 32,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -136,7 +146,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 33,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -148,24 +158,24 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 34,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Set Dirichlet boundary conditions.\n",
-    "dirichletBnd = ug4.DirichletBoundary2dCPU1()\n",
+    "dirichletBnd = ugcore.DirichletBoundary2dCPU1()\n",
     "dirichletBnd.add(1.0, \"u\", \"TOP_SC\")\n",
     "dirichletBnd.add(0.0, \"u\", \"BOTTOM_SC\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 35,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Create the global spatial discretization object.\n",
-    "domainDisc = ug4.DomainDiscretization2dCPU1(approxSpace)\n",
+    "domainDisc = ugcore.DomainDiscretization2dCPU1(approxSpace)\n",
     "domainDisc.add(elemDisc[\"LIP\"])\n",
     "domainDisc.add(elemDisc[\"COR\"])\n",
     "domainDisc.add(dirichletBnd)"
@@ -181,26 +191,35 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 36,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "ILUT: please use 'set_ordering_algorithm(..)' in the future\n"
+     ]
+    }
+   ],
    "source": [
-    "A = ug4.AssembledLinearOperatorCPU1(domainDisc)\n",
-    "u = ug4.GridFunction2dCPU1(approxSpace)\n",
-    "b = ug4.GridFunction2dCPU1(approxSpace)\n",
+    "A = ugcore.AssembledLinearOperatorCPU1(domainDisc)\n",
+    "u = ugcore.GridFunction2dCPU1(approxSpace)\n",
+    "b = ugcore.GridFunction2dCPU1(approxSpace)\n",
     "\n",
-    "ug4.Interpolate(0.0, u, \"u\")\n",
+    "ugcore.Interpolate(0.0, u, \"u\")\n",
     "\n",
     "domainDisc.assemble_linear(A, b)\n",
     "domainDisc.adjust_solution(u)\n",
     "\n",
-    "#lsolver = ug4.LUCPU1()\n",
-    "lsolver =slu.SuperLUCPU1()\n",
+    "lsolver = ugcore.LUCPU1()\n",
+    "#lsolver =slu.SuperLUCPU1()\n",
     "\n",
     "lsolver.init(A, u)\n",
     "lsolver.apply(u, b)\n",
     "\n",
-    "ug4.WriteGridFunctionToVTK(u, \"vtk/SkinSteadyState.vtk\")"
+    "ugcore.WriteGridFunctionToVTK(u, \"vtk/SkinSteadyState.vtk\")"
    ]
   },
   {
@@ -213,7 +232,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [
     {
@@ -222,20 +241,21 @@
      "text": [
       "Surface area [um^2]:\n",
       "30.100000000000016\n",
+      "IntegrateSubsets for subsets=BOTTOM_SC\n",
       "Surface fluxes [kg/s]:\n",
-      "-0.0030822575671149172\n",
-      "0.0030822575666958397\n"
+      "-0.0030822575664210083\n",
+      "0.0030822575668487456\n"
      ]
     }
    ],
    "source": [
-    "areaG=ug4.Integral(1.0, u, \"BOTTOM_SC\")\n",
+    "areaG=ugcore.Integral(1.0, u, \"BOTTOM_SC\")\n",
     "print(\"Surface area [um^2]:\")\n",
     "print(areaG)\n",
     "\n",
     "surfaceFlux = {}\n",
-    "surfaceFlux[\"BOT\"] = K[\"LIP\"]*D[\"LIP\"]*ug4.IntegrateNormalGradientOnManifold(u, \"u\", \"BOTTOM_SC\", \"LIP\")\n",
-    "surfaceFlux[\"TOP\"] = K[\"LIP\"]*D[\"LIP\"]*ug4.IntegrateNormalGradientOnManifold(u, \"u\", \"TOP_SC\", \"LIP\")\n",
+    "surfaceFlux[\"BOT\"] = K[\"LIP\"]*D[\"LIP\"]*ugcore.IntegrateNormalGradientOnManifold(u, \"u\", \"BOTTOM_SC\", \"LIP\")\n",
+    "surfaceFlux[\"TOP\"] = K[\"LIP\"]*D[\"LIP\"]*ugcore.IntegrateNormalGradientOnManifold(u, \"u\", \"TOP_SC\", \"LIP\")\n",
     "\n",
     "print(\"Surface fluxes [kg/s]:\")\n",
     "print(surfaceFlux[\"TOP\"])\n",
@@ -244,9 +264,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 38,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Normalized Fluxes [kg / (mu^2 * s)]:\n",
+      "-0.00010240058360202681\n",
+      "0.00010240058361623734\n"
+     ]
+    }
+   ],
    "source": [
     "print(\"Normalized Fluxes [kg / (mu^2 * s)]:\")\n",
     "print(surfaceFlux[\"TOP\"]/areaG)\n",
@@ -255,9 +285,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 39,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Relative Fluxes [1]:\n",
+      "-0.0018022502713956718\n",
+      "0.001802250271645777\n"
+     ]
+    }
+   ],
    "source": [
     "Jref = 1.0/17.6\n",
     "\n",
@@ -282,17 +322,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Create time discretization.\n",
-    "timeDisc=ug4.ThetaTimeStepCPU1(domainDisc, 1.0) "
+    "timeDisc=ugcore.ThetaTimeStepCPU1(domainDisc, 1.0) "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 41,
    "metadata": {
     "scrolled": true
    },
@@ -307,7 +347,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 42,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -337,7 +377,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 43,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -360,8 +400,8 @@
     "    \n",
     "    # 1) Compute fluxes.\n",
     "    gradFlux={}\n",
-    "    gradFlux[\"BOT\"] = ug4.IntegrateNormalGradientOnManifold(u, \"u\", \"BOTTOM_SC\", \"LIP\")\n",
-    "    gradFlux[\"TOP\"] = ug4.IntegrateNormalGradientOnManifold(u, \"u\", \"TOP_SC\", \"LIP\")\n",
+    "    gradFlux[\"BOT\"] = ugcore.IntegrateNormalGradientOnManifold(u, \"u\", \"BOTTOM_SC\", \"LIP\")\n",
+    "    gradFlux[\"TOP\"] = ugcore.IntegrateNormalGradientOnManifold(u, \"u\", \"TOP_SC\", \"LIP\")\n",
     "  \n",
     "    jTOP = K[\"LIP\"]*D[\"LIP\"]*gradFlux[\"TOP\"]\n",
     "    jBOT = K[\"LIP\"]*D[\"LIP\"]*gradFlux[\"BOT\"]\n",
@@ -384,12 +424,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 44,
    "metadata": {},
    "outputs": [],
    "source": [
     "# We use this function to define a callback\n",
-    "pyobserver= ug4.PythonCallbackObserver2dCPU1(MyPostProcess)"
+    "pyobserver= ugcore.PythonCallbackObserver2dCPU1(MyPostProcess)"
    ]
   },
   {
@@ -401,20 +441,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 45,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Create callback observer for file I/O.\n",
     "def MyVTKCallback(usol, step, time, dt) :\n",
-    "    ug4.WriteGridFunctionToVTK(usol, \"vtk/SkinDiffusion_\"+str(int(step)).zfill(5)+\".vtu\")\n",
+    "    ugcore.WriteGridFunctionToVTK(usol, \"vtk/SkinDiffusion_\"+str(int(step)).zfill(5)+\".vtu\")\n",
     "    \n",
-    "vtkobserver = ug4.PythonCallbackObserver2dCPU1(MyVTKCallback) "
+    "vtkobserver = ugcore.PythonCallbackObserver2dCPU1(MyVTKCallback) "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 46,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -432,7 +472,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 47,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -447,7 +487,7 @@
     "def MyProgressCallback(u, step, time, dt):\n",
     "    wProgress.value=time\n",
     "    \n",
-    "progressobserver = ug4.PythonCallbackObserver2dCPU1(MyProgressCallback)"
+    "progressobserver = ugcore.PythonCallbackObserver2dCPU1(MyProgressCallback)"
    ]
   },
   {
@@ -459,7 +499,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 48,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -478,19 +518,34 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 49,
    "metadata": {},
    "outputs": [],
    "source": [
     "# Set initial values.\n",
-    "ug4.Interpolate(0.0, u, \"u\")"
+    "ugcore.Interpolate(0.0, u, \"u\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 50,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "d006945a5bd14790a50691b81323970b",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "FloatProgress(value=0.0, description='Progress:', max=5000.0, style=ProgressStyle(bar_color='blue'))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
     "# Show progress bar.\n",
     "display(wProgress)"
@@ -498,9 +553,1020 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 51,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+++ Integrating: [\t0\t, \t5000\t] with dt=\t25(200 iters)\n",
+      "+++ Const timestep +++1(t=0, dt=25)\n",
+      "+++ Assemble (t=0, dt=25)\n",
+      "\n",
+      "ILUT: please use 'set_ordering_algorithm(..)' in the future\n",
+      "flux_top (\t25.0\t)=\t-0.10430773903686517\n",
+      "flux_bot (\t25.0\t)=\t3.2002967862279643e-11\n",
+      "mass_bot (\t25.0\t)=\t4.0003709827849554e-10\n",
+      "tlag=12.5\n",
+      "+++ Const timestep +++2(t=25, dt=25)\n",
+      "flux_top (\t50.0\t)=\t-0.033783767504077515\n",
+      "flux_bot (\t50.0\t)=\t3.012156241246096e-10\n",
+      "mass_bot (\t50.0\t)=\t4.565269498114611e-09\n",
+      "tlag=34.84384894913018\n",
+      "+++ Const timestep +++3(t=50, dt=25)\n",
+      "flux_top (\t75.0\t)=\t-0.027023337678278664\n",
+      "flux_bot (\t75.0\t)=\t1.5354299759488278e-09\n",
+      "mass_bot (\t75.0\t)=\t2.7523339499032578e-08\n",
+      "tlag=57.07450686116482\n",
+      "+++ Const timestep +++4(t=75, dt=25)\n",
+      "flux_top (\t100.0\t)=\t-0.023147365061062193\n",
+      "flux_bot (\t100.0\t)=\t5.6193937611384216e-09\n",
+      "mass_bot (\t100.0\t)=\t1.169586362126232e-07\n",
+      "tlag=79.18660959097325\n",
+      "+++ Const timestep +++5(t=100, dt=25)\n",
+      "flux_top (\t125.0\t)=\t-0.02056134363353459\n",
+      "flux_bot (\t125.0\t)=\t1.6530278292811834e-08\n",
+      "mass_bot (\t125.0\t)=\t3.938295368870014e-07\n",
+      "tlag=101.17526275657092\n",
+      "+++ Const timestep +++6(t=125, dt=25)\n",
+      "flux_top (\t150.0\t)=\t-0.01867607903956207\n",
+      "flux_bot (\t150.0\t)=\t4.1514067990049876e-08\n",
+      "mass_bot (\t150.0\t)=\t1.1193838654227728e-06\n",
+      "tlag=123.03603526180407\n",
+      "+++ Const timestep +++7(t=150, dt=25)\n",
+      "flux_top (\t175.0\t)=\t-0.017221420032417147\n",
+      "flux_bot (\t175.0\t)=\t9.237685929091508e-08\n",
+      "mass_bot (\t175.0\t)=\t2.7930204564348344e-06\n",
+      "tlag=144.7649337953892\n",
+      "+++ Const timestep +++8(t=175, dt=25)\n",
+      "flux_top (\t200.0\t)=\t-0.016054844059422865\n",
+      "flux_bot (\t200.0\t)=\t1.8672785522754273e-07\n",
+      "mass_bot (\t200.0\t)=\t6.281829387915557e-06\n",
+      "tlag=166.35837015178777\n",
+      "+++ Const timestep +++9(t=200, dt=25)\n",
+      "flux_top (\t225.0\t)=\t-0.01509286546727557\n",
+      "flux_bot (\t225.0\t)=\t3.490084775333214e-07\n",
+      "mass_bot (\t225.0\t)=\t1.297853354742636e-05\n",
+      "tlag=187.8131280960439\n",
+      "+++ Const timestep +++10(t=225, dt=25)\n",
+      "flux_top (\t250.0\t)=\t-0.01428257188245676\n",
+      "flux_bot (\t250.0\t)=\t6.11170191695567e-07\n",
+      "mass_bot (\t250.0\t)=\t2.4980766912787465e-05\n",
+      "tlag=209.1263329720263\n",
+      "+++ Const timestep +++11(t=250, dt=25)\n",
+      "flux_top (\t275.0\t)=\t-0.013588406530391967\n",
+      "flux_bot (\t275.0\t)=\t1.012911812612868e-06\n",
+      "mass_bot (\t275.0\t)=\t4.52817919666429e-05\n",
+      "tlag=230.2954251270545\n",
+      "+++ Const timestep +++12(t=275, dt=25)\n",
+      "flux_top (\t300.0\t)=\t-0.012985409862839704\n",
+      "flux_bot (\t300.0\t)=\t1.601442960262884e-06\n",
+      "mass_bot (\t300.0\t)=\t7.79612266275898e-05\n",
+      "tlag=251.31813710380786\n",
+      "+++ Const timestep +++13(t=300, dt=25)\n",
+      "flux_top (\t325.0\t)=\t-0.01245546669755071\n",
+      "flux_bot (\t325.0\t)=\t2.4307932909010255e-06\n",
+      "mass_bot (\t325.0\t)=\t0.00012836417976713868\n",
+      "tlag=272.1924740587227\n",
+      "+++ Const timestep +++14(t=325, dt=25)\n",
+      "flux_top (\t350.0\t)=\t-0.01198507620028769\n",
+      "flux_bot (\t350.0\t)=\t3.560729785602811e-06\n",
+      "mass_bot (\t350.0\t)=\t0.00020325821822343663\n",
+      "tlag=292.9166967273743\n",
+      "+++ Const timestep +++15(t=350, dt=25)\n",
+      "flux_top (\t375.0\t)=\t-0.011563953907444469\n",
+      "flux_bot (\t375.0\t)=\t5.055373058698394e-06\n",
+      "mass_bot (\t375.0\t)=\t0.0003109595037772017\n",
+      "tlag=313.48930629517884\n",
+      "+++ Const timestep +++16(t=375, dt=25)\n",
+      "flux_top (\t400.0\t)=\t-0.01118411810092349\n",
+      "flux_bot (\t400.0\t)=\t6.981617826106886e-06\n",
+      "mass_bot (\t400.0\t)=\t0.00046142188983726767\n",
+      "tlag=333.90903063873844\n",
+      "+++ Const timestep +++17(t=400, dt=25)\n",
+      "flux_top (\t425.0\t)=\t-0.010839272854492377\n",
+      "flux_bot (\t425.0\t)=\t9.407464031357028e-06\n",
+      "mass_bot (\t425.0\t)=\t0.0006662854130555666\n",
+      "tlag=354.17481152894135\n",
+      "+++ Const timestep +++18(t=425, dt=25)\n",
+      "flux_top (\t450.0\t)=\t-0.010524380240068333\n",
+      "flux_bot (\t450.0\t)=\t1.2400356572719076e-05\n",
+      "mass_bot (\t450.0\t)=\t0.0009388831706065179\n",
+      "tlag=374.2857925011551\n",
+      "+++ Const timestep +++19(t=450, dt=25)\n",
+      "flux_top (\t475.0\t)=\t-0.01023535704699207\n",
+      "flux_bot (\t475.0\t)=\t1.6025616412931643e-05\n",
+      "mass_bot (\t475.0\t)=\t0.0012942078329271519\n",
+      "tlag=394.241307193475\n",
+      "+++ Const timestep +++20(t=475, dt=25)\n",
+      "flux_top (\t500.0\t)=\t-0.009968855609882053\n",
+      "flux_bot (\t500.0\t)=\t2.034502720053598e-05\n",
+      "mass_bot (\t500.0\t)=\t0.0017488408780954971\n",
+      "tlag=414.04086802845734\n",
+      "+++ Const timestep +++21(t=500, dt=25)\n",
+      "flux_top (\t525.0\t)=\t-0.009722102705352696\n",
+      "flux_bot (\t525.0\t)=\t2.5415621978224123e-05\n",
+      "mass_bot (\t525.0\t)=\t0.0023208489928299983\n",
+      "tlag=433.6841551696637\n",
+      "+++ Const timestep +++22(t=525, dt=25)\n",
+      "flux_top (\t550.0\t)=\t-0.009492779303002878\n",
+      "flux_bot (\t550.0\t)=\t3.12886960195957e-05\n",
+      "mass_bot (\t550.0\t)=\t0.0030296529678027462\n",
+      "tlag=453.17100572343077\n",
+      "+++ Const timestep +++23(t=550, dt=25)\n",
+      "flux_top (\t575.0\t)=\t-0.009278929553457438\n",
+      "flux_bot (\t575.0\t)=\t3.800905557422246e-05\n",
+      "mass_bot (\t575.0\t)=\t0.0038958748627254732\n",
+      "tlag=472.50140318225544\n",
+      "+++ Const timestep +++24(t=575, dt=25)\n",
+      "flux_top (\t600.0\t)=\t-0.00907889101570265\n",
+      "flux_bot (\t600.0\t)=\t4.561449898301686e-05\n",
+      "mass_bot (\t600.0\t)=\t0.004941169294690965\n",
+      "tlag=491.6754671243752\n",
+      "+++ Const timestep +++25(t=600, dt=25)\n",
+      "flux_top (\t625.0\t)=\t-0.008891240537318724\n",
+      "flux_bot (\t625.0\t)=\t5.413551646210508e-05\n",
+      "mass_bot (\t625.0\t)=\t0.006188044487754989\n",
+      "tlag=510.6934431929429\n",
+      "+++ Const timestep +++26(t=625, dt=25)\n",
+      "flux_top (\t650.0\t)=\t-0.008714751817590282\n",
+      "flux_bot (\t650.0\t)=\t6.359518772820015e-05\n",
+      "mass_bot (\t650.0\t)=\t0.007659678290133805\n",
+      "tlag=529.5556933824844\n",
+      "+++ Const timestep +++27(t=650, dt=25)\n",
+      "flux_top (\t675.0\t)=\t-0.008548361796210328\n",
+      "flux_bot (\t675.0\t)=\t7.400925222770833e-05\n",
+      "mass_bot (\t675.0\t)=\t0.009379733789582661\n",
+      "tlag=548.262686660804\n",
+      "+++ Const timestep +++28(t=675, dt=25)\n",
+      "flux_top (\t700.0\t)=\t-0.008391143786298304\n",
+      "flux_bot (\t700.0\t)=\t8.538632460222803e-05\n",
+      "mass_bot (\t700.0\t)=\t0.011372178499956864\n",
+      "tlag=566.814989953788\n",
+      "+++ Const timestep +++29(t=700, dt=25)\n",
+      "flux_top (\t725.0\t)=\t-0.008242285809167213\n",
+      "flux_bot (\t725.0\t)=\t9.772822770091635e-05\n",
+      "mass_bot (\t725.0\t)=\t0.013661110403746169\n",
+      "tlag=585.2132595143944\n",
+      "+++ Const timestep +++30(t=725, dt=25)\n",
+      "flux_top (\t750.0\t)=\t-0.008101072989034714\n",
+      "flux_bot (\t750.0\t)=\t0.00011103041648384794\n",
+      "mass_bot (\t750.0\t)=\t0.01627059345605572\n",
+      "tlag=603.4582326958787\n",
+      "+++ Const timestep +++31(t=750, dt=25)\n",
+      "flux_top (\t775.0\t)=\t-0.007966873133943588\n",
+      "flux_bot (\t775.0\t)=\t0.00012528246813254884\n",
+      "mass_bot (\t775.0\t)=\t0.01922450451376068\n",
+      "tlag=621.550720142094\n",
+      "+++ Const timestep +++32(t=775, dt=25)\n",
+      "flux_top (\t800.0\t)=\t-0.007839124849690434\n",
+      "flux_bot (\t800.0\t)=\t0.00014046861624748402\n",
+      "mass_bot (\t800.0\t)=\t0.02254639306851109\n",
+      "tlag=639.4915984023945\n",
+      "+++ Const timestep +++33(t=800, dt=25)\n",
+      "flux_top (\t825.0\t)=\t-0.007717327673696276\n",
+      "flux_bot (\t825.0\t)=\t0.00015656830988676455\n",
+      "mass_bot (\t825.0\t)=\t0.0262593546451892\n",
+      "tlag=657.2818029767272\n",
+      "+++ Const timestep +++34(t=825, dt=25)\n",
+      "flux_top (\t850.0\t)=\t-0.0076010338403198485\n",
+      "flux_bot (\t850.0\t)=\t0.00017355678115956198\n",
+      "mass_bot (\t850.0\t)=\t0.03038591828326828\n",
+      "tlag=674.9223217885532\n",
+      "+++ Const timestep +++35(t=850, dt=25)\n",
+      "flux_top (\t875.0\t)=\t-0.00748984136302667\n",
+      "flux_bot (\t875.0\t)=\t0.00019140560798035355\n",
+      "mass_bot (\t875.0\t)=\t0.03494794814751723\n",
+      "tlag=692.4141890811037\n",
+      "flux_top (\t900.0\t)=\t-0.007383388195876656\n",
+      "flux_bot (\t900.0\t)=\t0.00021008326130294831\n",
+      "mass_bot (\t900.0\t)=\t0.0399665590135585\n",
+      "tlag=709.7584797299717\n",
+      "+++ Const timestep +++36(t=875, dt=25)\n",
+      "+++ Const timestep +++37(t=900, dt=25)\n",
+      "flux_top (\t925.0\t)=\t-0.007281347282051364\n",
+      "flux_bot (\t925.0\t)=\t0.00022955562861390962\n",
+      "mass_bot (\t925.0\t)=\t0.045462045137519225\n",
+      "tlag=726.9563039598477\n",
+      "flux_top (\t950.0\t)=\t-0.007183422328593999\n",
+      "flux_bot (\t950.0\t)=\t0.0002497865076426863\n",
+      "mass_bot (\t950.0\t)=\t0.05145382184072667\n",
+      "tlag=744.0088024516915\n",
+      "+++ Const timestep +++38(t=925, dt=25)\n",
+      "+++ Const timestep +++39(t=950, dt=25)\n",
+      "flux_top (\t975.0\t)=\t-0.007089344191786943\n",
+      "flux_bot (\t975.0\t)=\t0.0002707380661269617\n",
+      "mass_bot (\t975.0\t)=\t0.05796037901284727\n",
+      "tlag=760.9171418264215\n",
+      "+++ Const timestep +++40(t=975, dt=25)\n",
+      "flux_top (\t1000.0\t)=\t-0.006998867766962409\n",
+      "flux_bot (\t1000.0\t)=\t0.00029237126505360426\n",
+      "mass_bot (\t1000.0\t)=\t0.06499924565260434\n",
+      "tlag=777.6825104864968\n",
+      "+++ Const timestep +++41(t=1000, dt=25)\n",
+      "flux_top (\t1025.0\t)=\t-0.006911769304446682\n",
+      "flux_bot (\t1025.0\t)=\t0.0003146462441087541\n",
+      "mass_bot (\t1025.0\t)=\t0.07258696451713383\n",
+      "tlag=794.3061147996893\n",
+      "+++ Const timestep +++42(t=1025, dt=25)\n",
+      "flux_top (\t1050.0\t)=\t-0.006827844082454649\n",
+      "flux_bot (\t1050.0\t)=\t0.000337522669121494\n",
+      "mass_bot (\t1050.0\t)=\t0.08073907593251192\n",
+      "tlag=810.7891756051229\n",
+      "flux_top (\t1075.0\t)=\t-0.00674690438258159+++ Const timestep +++43(t=1050, dt=25)\n",
+      "\n",
+      "flux_bot (\t1075.0\t)=\t0.000360960042116986\n",
+      "mass_bot (\t1075.0\t)=\t0.08947010982299292\n",
+      "tlag=827.13292502333\n",
+      "+++ Const timestep +++44(t=1075, dt=25)\n",
+      "flux_top (\t1100.0\t)=\t-0.006668777724165765\n",
+      "flux_bot (\t1100.0\t)=\t0.00038491797522730497\n",
+      "mass_bot (\t1100.0\t)=\t0.09879358503979656\n",
+      "tlag=843.3386035519486\n",
+      "+++ Const timestep +++45(t=1100, dt=25)\n",
+      "flux_top (\t1125.0\t)=\t-0.006593305316672153\n",
+      "flux_bot (\t1125.0\t)=\t0.0004093564301718573\n",
+      "mass_bot (\t1125.0\t)=\t0.10872201510728609\n",
+      "tlag=859.4074574286226\n",
+      "+++ Const timestep +++46(t=1125, dt=25)\n",
+      "flux_top (\t1150.0\t)=\t-0.006520340700457898\n",
+      "flux_bot (\t1150.0\t)=\t0.00043423592534118014\n",
+      "mass_bot (\t1150.0\t)=\t0.11926691955119906\n",
+      "tlag=875.3407362426525\n",
+      "+++ Const timestep +++47(t=1150, dt=25)\n",
+      "flux_top (\t1175.0\t)=\t-0.006449748550277965\n",
+      "flux_bot (\t1175.0\t)=\t0.0004595177127238968\n",
+      "mass_bot (\t1175.0\t)=\t0.13043884002701253\n",
+      "tlag=891.139690777519\n",
+      "+++ Const timestep +++48(t=1175, dt=25)\n",
+      "flux_top (\t1200.0\t)=\t-0.006381403615853403\n",
+      "flux_bot (\t1200.0\t)=\t0.00048516392703089995\n",
+      "mass_bot (\t1200.0\t)=\t0.14224736052394749\n",
+      "tlag=906.8055710685849\n",
+      "+++ Const timestep +++49(t=1200, dt=25)\n",
+      "flux_top (\t1225.0\t)=\t-0.006315189785961283\n",
+      "flux_bot (\t1225.0\t)=\t0.0005111377094033086\n",
+      "mass_bot (\t1225.0\t)=\t0.1547011309793751\n",
+      "tlag=922.3396246581574\n",
+      "+++ Const timestep +++50(t=1225, dt=25)\n",
+      "flux_top (\t1250.0\t)=\t-0.006250999254356142\n",
+      "flux_bot (\t1250.0\t)=\t0.0005374033080755061\n",
+      "mass_bot (\t1250.0\t)=\t0.16780789369786028\n",
+      "tlag=937.7430950345341\n",
+      "+++ Const timestep +++51(t=1250, dt=25)\n",
+      "flux_top (\t1275.0\t)=\t-0.006188731777938728\n",
+      "flux_bot (\t1275.0\t)=\t0.0005639261582915338\n",
+      "mass_bot (\t1275.0\t)=\t0.1815745120274483\n",
+      "tlag=953.017220237584\n",
+      "+++ Const timestep +++52(t=1275, dt=25)\n",
+      "flux_top (\t1300.0\t)=\t-0.006128294015723932\n",
+      "flux_bot (\t1300.0\t)=\t0.0005906729436860929\n",
+      "mass_bot (\t1300.0\t)=\t0.19600700080216862\n",
+      "tlag=968.1632316202135\n",
+      "+++ Const timestep +++53(t=1300, dt=25)\n",
+      "flux_top (\t1325.0\t)=\t-0.006069598934574213\n",
+      "flux_bot (\t1325.0\t)=\t0.0006176116412132823\n",
+      "mass_bot (\t1325.0\t)=\t0.2111105581134108\n",
+      "tlag=983.1823527505254\n",
+      "+++ Const timestep +++54(t=1325, dt=25)\n",
+      "flux_top (\t1350.0\t)=\t-0.006012565278341164\n",
+      "flux_bot (\t1350.0\t)=\t0.0006447115515790183\n",
+      "mass_bot (\t1350.0\t)=\t0.22688959802331457\n",
+      "tlag=998.0757984440952\n",
+      "+++ Const timestep +++55(t=1350, dt=25)\n",
+      "flux_top (\t1375.0\t)=\t-0.005957117089951653\n",
+      "flux_bot (\t1375.0\t)=\t0.000671943316988093\n",
+      "mass_bot (\t1375.0\t)=\t0.24334778388040346\n",
+      "tlag=1012.8447739145895\n",
+      "+++ Const timestep +++56(t=1375, dt=25)\n",
+      "flux_top (\t1400.0\t)=\t-0.005903183280601442\n",
+      "flux_bot (\t1400.0\t)=\t0.0006992789278695039\n",
+      "mass_bot (\t1400.0\t)=\t0.2604880619411234\n",
+      "tlag=1027.4904740305083\n",
+      "+++ Const timestep +++57(t=1400, dt=25)\n",
+      "flux_top (\t1425.0\t)=\t-0.0058506972418688775\n",
+      "flux_bot (\t1425.0\t)=\t0.0007266917201109892\n",
+      "mass_bot (\t1425.0\t)=\t0.2783126950408796\n",
+      "tlag=1042.014082672674\n",
+      "+++ Const timestep +++58(t=1425, dt=25)\n",
+      "flux_top (\t1450.0\t)=\t-0.005799596493627673\n",
+      "flux_bot (\t1450.0\t)=\t0.0007541563641769893\n",
+      "mass_bot (\t1450.0\t)=\t0.2968232960944793\n",
+      "tlag=1056.416772179066\n",
+      "+++ Const timestep +++59(t=1450, dt=25)\n",
+      "flux_top (\t1475.0\t)=\t-0.005749822366261271\n",
+      "flux_bot (\t1475.0\t)=\t0.0007816488473604733\n",
+      "mass_bot (\t1475.0\t)=\t0.3160208612386976\n",
+      "tlag=1070.6997028705932\n",
+      "+++ Const timestep +++60(t=1475, dt=25)\n",
+      "flux_top (\t1500.0\t)=\t-0.005701319710705217\n",
+      "flux_bot (\t1500.0\t)=\t0.0008091464502900772\n",
+      "mass_bot (\t1500.0\t)=\t0.3359058024593295\n",
+      "tlag=1084.8640226513903\n",
+      "+++ Const timestep +++61(t=1500, dt=25)\n",
+      "flux_top (\t1525.0\t)=\t-0.005654036636980084\n",
+      "flux_bot (\t1525.0\t)=\t0.0008366277186861512\n",
+      "mass_bot (\t1525.0\t)=\t0.3564779795715324\n",
+      "tlag=1098.9108666739503\n",
+      "+++ Const timestep +++62(t=1525, dt=25)\n",
+      "flux_top (\t1550.0\t)=\t-0.005607924273807367\n",
+      "flux_bot (\t1550.0\t)=\t0.0008640724312574539\n",
+      "mass_bot (\t1550.0\t)=\t0.3777367314458274\n",
+      "tlag=1112.8413570652629\n",
+      "+++ Const timestep +++63(t=1550, dt=25)\n",
+      "flux_top (\t1575.0\t)=\t-0.005562936550406367\n",
+      "flux_bot (\t1575.0\t)=\t0.0008914615645185932\n",
+      "mass_bot (\t1575.0\t)=\t0.399680906393028\n",
+      "tlag=1126.6566027062945\n",
+      "+++ Const timestep +++64(t=1575, dt=25)\n",
+      "flux_top (\t1600.0\t)=\t-0.005519029997278649\n",
+      "flux_bot (\t1600.0\t)=\t0.0009187772552202232\n",
+      "mass_bot (\t1600.0\t)=\t0.4223088916397632\n",
+      "tlag=1140.357699061086\n",
+      "+++ Const timestep +++65(t=1600, dt=25)\n",
+      "flux_top (\t1625.0\t)=\t-0.005476163563023346\n",
+      "flux_bot (\t1625.0\t)=\t0.0009460027609891695\n",
+      "mass_bot (\t1625.0\t)=\t0.4456186418423806\n",
+      "tlag=1153.945728047952\n",
+      "+++ Const timestep +++66(t=1625, dt=25)\n",
+      "flux_top (\t1650.0\t)=\t-0.005434298447500318\n",
+      "flux_bot (\t1650.0\t)=\t0.0009731224197088896\n",
+      "mass_bot (\t1650.0\t)=\t0.46960770660110635\n",
+      "tlag=1167.421757951492\n",
+      "+++ Const timestep +++67(t=1650, dt=25)\n",
+      "flux_top (\t1675.0\t)=\t-0.005393397947327992\n",
+      "flux_bot (\t1675.0\t)=\t0.0010001216080894862\n",
+      "mass_bot (\t1675.0\t)=\t0.49427325694858604\n",
+      "tlag=1180.7868433691908\n",
+      "+++ Const timestep +++68(t=1675, dt=25)\n",
+      "flux_top (\t1700.0\t)=\t-0.005353427315043953\n",
+      "flux_bot (\t1700.0\t)=\t0.0010269866998182113\n",
+      "mass_bot (\t1700.0\t)=\t0.5196121107974323\n",
+      "tlag=1194.042025189411\n",
+      "+++ Const timestep +++69(t=1700, dt=25)\n",
+      "flux_top (\t1725.0\t)=\t-0.005314353629956969\n",
+      "flux_bot (\t1725.0\t)=\t0.0010537050236249824\n",
+      "mass_bot (\t1725.0\t)=\t0.5456207573404722\n",
+      "tlag=1207.188330598051\n",
+      "+++ Const timestep +++70(t=1725, dt=25)\n",
+      "flux_top (\t1750.0\t)=\t-0.005276145676594102\n",
+      "flux_bot (\t1750.0\t)=\t0.0010802648215428058\n",
+      "mass_bot (\t1750.0\t)=\t0.5722953804050696\n",
+      "tlag=1220.2267731094562\n",
+      "+++ Const timestep +++71(t=1750, dt=25)\n",
+      "flux_top (\t1775.0\t)=\t-0.0052387738356784045\n",
+      "flux_bot (\t1775.0\t)=\t0.0011066552076030424\n",
+      "mass_bot (\t1775.0\t)=\t0.5996318807693927\n",
+      "tlag=1233.1583526199058\n",
+      "+++ Const timestep +++72(t=1775, dt=25)\n",
+      "flux_top (\t1800.0\t)=\t-0.005202209981335462\n",
+      "flux_bot (\t1800.0\t)=\t0.0011328661271635364\n",
+      "mass_bot (\t1800.0\t)=\t0.6276258974539749\n",
+      "tlag=1245.9840554810999\n",
+      "+++ Const timestep +++73(t=1800, dt=25)\n",
+      "flux_top (\t1825.0\t)=\t-0.005166427385907689\n",
+      "flux_bot (\t1825.0\t)=\t0.00115888831703257\n",
+      "mass_bot (\t1825.0\t)=\t0.6562728280064263\n",
+      "tlag=1258.7048545912799\n",
+      "+++ Const timestep +++74(t=1825, dt=25)\n",
+      "flux_top (\t1850.0\t)=\t-0.0051314006344912065\n",
+      "flux_bot (\t1850.0\t)=\t0.0011847132665188597\n",
+      "mass_bot (\t1850.0\t)=\t0.6855678478008191\n",
+      "tlag=1271.3217095008317\n",
+      "+++ Const timestep +++75(t=1850, dt=25)\n",
+      "flux_top (\t1875.0\t)=\t-0.005097105541986042\n",
+      "flux_bot (\t1875.0\t)=\t0.0012103331795159493\n",
+      "mass_bot (\t1875.0\t)=\t0.7155059283762543\n",
+      "tlag=1283.8355665318472\n",
+      "+++ Const timestep +++76(t=1875, dt=25)\n",
+      "flux_top (\t1900.0\t)=\t-0.005063519078691691\n",
+      "flux_bot (\t1900.0\t)=\t0.001235740937705849\n",
+      "mass_bot (\t1900.0\t)=\t0.7460818548415268\n",
+      "tlag=1296.2473589111432\n",
+      "+++ Const timestep +++77(t=1900, dt=25)\n",
+      "flux_top (\t1925.0\t)=\t-0.00503061930048655\n",
+      "flux_bot (\t1925.0\t)=\t0.0012609300649327603\n",
+      "mass_bot (\t1925.0\t)=\t0.7772902423745094\n",
+      "tlag=1308.5580069097973\n",
+      "+++ Const timestep +++78(t=1925, dt=25)\n",
+      "flux_top (\t1950.0\t)=\t-0.004998385283462745\n",
+      "flux_bot (\t1950.0\t)=\t0.0012858946928080636\n",
+      "mass_bot (\t1950.0\t)=\t0.8091255518462697\n",
+      "tlag=1320.768417995919\n",
+      "+++ Const timestep +++79(t=1950, dt=25)\n",
+      "flux_top (\t1975.0\t)=\t-0.004966797064305438\n",
+      "flux_bot (\t1975.0\t)=\t0.0013106295275591812\n",
+      "mass_bot (\t1975.0\t)=\t0.8415821046008602\n",
+      "tlag=1332.8794869911408\n",
+      "+++ Const timestep +++80(t=1975, dt=25)\n",
+      "flux_top (\t2000.0\t)=\t-0.0049358355831493434\n",
+      "flux_bot (\t2000.0\t)=\t0.0013351298181559325\n",
+      "mass_bot (\t2000.0\t)=\t0.8746540964222992\n",
+      "tlag=1344.8920962379807\n",
+      "+++ Const timestep +++81(t=2000, dt=25)\n",
+      "flux_top (\t2025.0\t)=\t-0.004905482632051217\n",
+      "flux_bot (\t2025.0\t)=\t0.0013593913257029027\n",
+      "mass_bot (\t2025.0\t)=\t0.9083356107205346\n",
+      "tlag=1356.8071157686254\n",
+      "+++ Const timestep +++82(t=2025, dt=25)\n",
+      "flux_top (\t2050.0\t)=\t-0.004875720805388788\n",
+      "flux_bot (\t2050.0\t)=\t0.0013834102941051734\n",
+      "mass_bot (\t2050.0\t)=\t0.9426206309681355\n",
+      "tlag=1368.6254034795602\n",
+      "+++ Const timestep +++83(t=2050, dt=25)\n",
+      "flux_top (\t2075.0\t)=\t-0.004846533454070125\n",
+      "flux_bot (\t2075.0\t)=\t0.0014071834220004146\n",
+      "mass_bot (\t2075.0\t)=\t0.9775030524194553\n",
+      "tlag=1380.34780531321\n",
+      "+++ Const timestep +++84(t=2075, dt=25)\n",
+      "flux_top (\t2100.0\t)=\t-0.004817904643905116\n",
+      "flux_bot (\t2100.0\t)=\t0.0014307078359237682\n",
+      "mass_bot (\t2100.0\t)=\t1.0129766931435076\n",
+      "tlag=1391.9751554380377\n",
+      "+++ Const timestep +++85(t=2100, dt=25)\n",
+      "flux_top (\t2125.0\t)=\t-0.004789819114998177\n",
+      "flux_bot (\t2125.0\t)=\t0.0014539810646977424\n",
+      "mass_bot (\t2125.0\t)=\t1.0490353044012766\n",
+      "tlag=1403.50827643388\n",
+      "+++ Const timestep +++86(t=2125, dt=25)\n",
+      "flux_top (\t2150.0\t)=\t-0.004762262244211429\n",
+      "flux_bot (\t2150.0\t)=\t0.0014770010150120944\n",
+      "mass_bot (\t2150.0\t)=\t1.0856725803976495\n",
+      "tlag=1414.947979478024\n",
+      "+++ Const timestep +++87(t=2150, dt=25)\n",
+      "flux_top (\t2175.0\t)=\t-0.004735220011223162\n",
+      "flux_bot (\t2175.0\t)=\t0.00149976594816531\n",
+      "mass_bot (\t2175.0\t)=\t1.122882167437367\n",
+      "tlag=1426.2950645325636\n",
+      "+++ Const timestep +++88(t=2175, dt=25)\n",
+      "flux_top (\t2200.0\t)=\t-0.004708678964437349\n",
+      "flux_bot (\t2200.0\t)=\t0.0015222744579390247\n",
+      "mass_bot (\t2200.0\t)=\t1.1606576725136712\n",
+      "tlag=1437.5503205347995\n",
+      "+++ Const timestep +++89(t=2200, dt=25)\n",
+      "flux_top (\t2225.0\t)=\t-0.00468262619220275\n",
+      "flux_bot (\t2225.0\t)=\t0.0015445254495586665\n",
+      "mass_bot (\t2225.0\t)=\t1.1989926713573924\n",
+      "tlag=1448.7145255845458\n",
+      "+++ Const timestep +++90(t=2225, dt=25)\n",
+      "flux_top (\t2250.0\t)=\t-0.004657049292450608\n",
+      "flux_bot (\t2250.0\t)=\t0.0015665181197170594\n",
+      "mass_bot (\t2250.0\t)=\t1.237880715973339\n",
+      "tlag=1459.7884471346415\n",
+      "+++ Const timestep +++91(t=2250, dt=25)\n",
+      "flux_top (\t2275.0\t)=\t-0.004631936346760955\n",
+      "flux_bot (\t2275.0\t)=\t0.0015882519376136392\n",
+      "mass_bot (\t2275.0\t)=\t1.2773153416899727\n",
+      "tlag=1470.7728421794664\n",
+      "+++ Const timestep +++92(t=2275, dt=25)\n",
+      "flux_top (\t2300.0\t)=\t-0.004607275894141097\n",
+      "flux_bot (\t2300.0\t)=\t0.0016097266269734002\n",
+      "mass_bot (\t2300.0\t)=\t1.3172900737473108\n",
+      "tlag=1481.6684574423218\n",
+      "+++ Const timestep +++93(t=2300, dt=25)\n",
+      "flux_top (\t2325.0\t)=\t-0.0045830569082724254\n",
+      "flux_bot (\t2325.0\t)=\t0.0016309421490134735\n",
+      "mass_bot (\t2325.0\t)=\t1.3577984334471467\n",
+      "tlag=1492.4760295645965\n",
+      "+++ Const timestep +++94(t=2325, dt=25)\n",
+      "flux_top (\t2350.0\t)=\t-0.004559268774604802\n",
+      "flux_bot (\t2350.0\t)=\t0.0016518986863087396\n",
+      "mass_bot (\t2350.0\t)=\t1.3988339438886743\n",
+      "tlag=1503.1962852912927\n",
+      "+++ Const timestep +++95(t=2350, dt=25)\n",
+      "flux_top (\t2375.0\t)=\t-0.004535901269463083\n",
+      "flux_bot (\t2375.0\t)=\t0.0016725966275261768\n",
+      "mass_bot (\t2375.0\t)=\t1.4403901353116109\n",
+      "tlag=1513.8299416566483\n",
+      "+++ Const timestep +++96(t=2375, dt=25)\n",
+      "flux_top (\t2400.0\t)=\t-0.004512944540295173\n",
+      "flux_bot (\t2400.0\t)=\t0.0016930365529857763\n",
+      "mass_bot (\t2400.0\t)=\t1.4824605500680104\n",
+      "tlag=1524.3777061673015\n",
+      "+++ Const timestep +++97(t=2400, dt=25)\n",
+      "flux_top (\t2425.0\t)=\t-0.004490389087598633\n",
+      "flux_bot (\t2425.0\t)=\t0.001713219221016791\n",
+      "mass_bot (\t2425.0\t)=\t1.5250387472430424\n",
+      "tlag=1534.8402769857228\n",
+      "+++ Const timestep +++98(t=2425, dt=25)\n",
+      "flux_top (\t2450.0\t)=\t-0.004468225745986305\n",
+      "flux_bot (\t2450.0\t)=\t0.0017331455550637688\n",
+      "mass_bot (\t2450.0\t)=\t1.5681183069440494\n",
+      "tlag=1545.2183431088956\n",
+      "+++ Const timestep +++99(t=2450, dt=25)\n",
+      "flux_top (\t2475.0\t)=\t-0.0044464456703143165\n",
+      "flux_bot (\t2475.0\t)=\t0.0017528166315187582\n",
+      "mass_bot (\t2475.0\t)=\t1.611692834276331\n",
+      "tlag=1555.5125845480757\n",
+      "+++ Const timestep +++100(t=2475, dt=25)\n",
+      "flux_top (\t2500.0\t)=\t-0.0044250403167686875\n",
+      "flux_bot (\t2500.0\t)=\t0.0017722336682366043\n",
+      "mass_bot (\t2500.0\t)=\t1.655755963023273\n",
+      "tlag=1565.7236725048951\n",
+      "+++ Const timestep +++101(t=2500, dt=25)\n",
+      "flux_top (\t2525.0\t)=\t-0.004404001431970831\n",
+      "flux_bot (\t2525.0\t)=\t0.0017913980137056014\n",
+      "mass_bot (\t2525.0\t)=\t1.7003013590475506\n",
+      "tlag=1575.852269546516\n",
+      "+++ Const timestep +++102(t=2525, dt=25)\n",
+      "flux_top (\t2550.0\t)=\t-0.004383321035285378\n",
+      "flux_bot (\t2550.0\t)=\t0.0018103111368356999\n",
+      "mass_bot (\t2550.0\t)=\t1.745322723429317\n",
+      "tlag=1585.8990297767145\n",
+      "+++ Const timestep +++103(t=2550, dt=25)\n",
+      "flux_top (\t2575.0\t)=\t-0.004362991409736293\n",
+      "flux_bot (\t2575.0\t)=\t0.0018289746173399556\n",
+      "mass_bot (\t2575.0\t)=\t1.7908137953565126\n",
+      "tlag=1595.8645990062691\n",
+      "+++ Const timestep +++104(t=2575, dt=25)\n",
+      "flux_top (\t2600.0\t)=\t-0.004343005086529274\n",
+      "flux_bot (\t2600.0\t)=\t0.0018473901366778663\n",
+      "mass_bot (\t2600.0\t)=\t1.8367683547817355\n",
+      "tlag=1605.7496149217468\n",
+      "+++ Const timestep +++105(t=2600, dt=25)\n",
+      "flux_top (\t2625.0\t)=\t-0.004323354835145942\n",
+      "flux_bot (\t2625.0\t)=\t0.0018655594695267035\n",
+      "mass_bot (\t2625.0\t)=\t1.8831802248592926\n",
+      "tlag=1615.554707249799\n",
+      "+++ Const timestep +++106(t=2625, dt=25)\n",
+      "flux_top (\t2650.0\t)=\t-0.0043040336528336704\n",
+      "flux_bot (\t2650.0\t)=\t0.0018834844757574802\n",
+      "mass_bot (\t2650.0\t)=\t1.9300432741753448\n",
+      "tlag=1625.2804979190817\n",
+      "+++ Const timestep +++107(t=2650, dt=25)\n",
+      "flux_top (\t2675.0\t)=\t-0.004285034752349022\n",
+      "flux_bot (\t2675.0\t)=\t0.0019011670928959531\n",
+      "mass_bot (\t2675.0\t)=\t1.9773514187835128\n",
+      "tlag=1634.9276012233559\n",
+      "+++ Const timestep +++108(t=2675, dt=25)\n",
+      "flux_top (\t2700.0\t)=\t-0.00426635155395353\n",
+      "flux_bot (\t2700.0\t)=\t0.0019186093290300427\n",
+      "mass_bot (\t2700.0\t)=\t2.0250986240575877\n",
+      "tlag=1644.496623978483\n",
+      "+++ Const timestep +++109(t=2700, dt=25)\n",
+      "flux_top (\t2725.0\t)=\t-0.004247977675687339\n",
+      "flux_bot (\t2725.0\t)=\t0.0019358132561452695\n",
+      "mass_bot (\t2725.0\t)=\t2.073278906372279\n",
+      "tlag=1653.9881656763002\n",
+      "+++ Const timestep +++110(t=2725, dt=25)\n",
+      "flux_top (\t2750.0\t)=\t-0.004229906924611248\n",
+      "flux_bot (\t2750.0\t)=\t0.001952781003874651\n",
+      "mass_bot (\t2750.0\t)=\t2.1218863346225283\n",
+      "tlag=1663.4028186405217\n",
+      "+++ Const timestep +++111(t=2750, dt=25)\n",
+      "flux_top (\t2775.0\t)=\t-0.0042121332879845085\n",
+      "flux_bot (\t2775.0\t)=\t0.001969514753625682\n",
+      "mass_bot (\t2775.0\t)=\t2.1709150315912824\n",
+      "tlag=1672.7411681761498\n",
+      "+++ Const timestep +++112(t=2775, dt=25)\n",
+      "flux_top (\t2800.0\t)=\t-0.004194650926768352\n",
+      "flux_bot (\t2800.0\t)=\t0.0019860167330698494\n",
+      "mass_bot (\t2800.0\t)=\t2.2203591751749765\n",
+      "tlag=1682.0037927158364\n",
+      "+++ Const timestep +++113(t=2800, dt=25)\n",
+      "flux_top (\t2825.0\t)=\t-0.004177454167223071\n",
+      "flux_bot (\t2825.0\t)=\t0.0020022892109827165\n",
+      "mass_bot (\t2825.0\t)=\t2.2702129994756337\n",
+      "tlag=1691.191263967596\n",
+      "+++ Const timestep +++114(t=2825, dt=25)\n",
+      "flux_top (\t2850.0\t)=\t-0.004160537493582315\n",
+      "flux_bot (\t2850.0\t)=\t0.0020183344924028074\n",
+      "mass_bot (\t2850.0\t)=\t2.3204707957679527\n",
+      "tlag=1700.304147056688\n",
+      "+++ Const timestep +++115(t=2850, dt=25)\n",
+      "flux_top (\t2875.0\t)=\t-0.004143895543202833\n",
+      "flux_bot (\t2875.0\t)=\t0.00203415491409985\n",
+      "mass_bot (\t2875.0\t)=\t2.371126913349236\n",
+      "tlag=1709.343000666445\n",
+      "+++ Const timestep +++116(t=2875, dt=25)\n",
+      "flux_top (\t2900.0\t)=\t-0.004127523098253175\n",
+      "flux_bot (\t2900.0\t)=\t0.002049752840329371\n",
+      "mass_bot (\t2900.0\t)=\t2.4221757602796012\n",
+      "tlag=1718.3083771746908\n",
+      "+++ Const timestep +++117(t=2900, dt=25)\n",
+      "flux_top (\t2925.0\t)=\t-0.0041114150807203455\n",
+      "flux_bot (\t2925.0\t)=\t0.0020651306588665647\n",
+      "mass_bot (\t2925.0\t)=\t2.4736118040195505\n",
+      "tlag=1727.2008227909519\n",
+      "+++ Const timestep +++118(t=2925, dt=25)\n",
+      "flux_top (\t2950.0\t)=\t-0.004095566547095993\n",
+      "flux_bot (\t2950.0\t)=\t0.0020802907772928835\n",
+      "mass_bot (\t2950.0\t)=\t2.525429571971544\n",
+      "tlag=1736.0208776881054\n",
+      "+++ Const timestep +++119(t=2950, dt=25)\n",
+      "flux_top (\t2975.0\t)=\t-0.004079972682578039\n",
+      "flux_bot (\t2975.0\t)=\t0.002095235619527811\n",
+      "mass_bot (\t2975.0\t)=\t2.5776236519318023\n",
+      "tlag=1744.7690761324955\n",
+      "+++ Const timestep +++120(t=2975, dt=25)\n",
+      "flux_top (\t3000.0\t)=\t-0.004064628796141208\n",
+      "flux_bot (\t3000.0\t)=\t0.0021099676225906144\n",
+      "mass_bot (\t3000.0\t)=\t2.6301886924582827\n",
+      "tlag=1753.445946611758\n",
+      "+++ Const timestep +++121(t=3000, dt=25)\n",
+      "flux_top (\t3025.0\t)=\t-0.004049530314807321\n",
+      "flux_bot (\t3025.0\t)=\t0.0021244892335818256\n",
+      "mass_bot (\t3025.0\t)=\t2.683119403160438\n",
+      "tlag=1762.0520119620571\n",
+      "+++ Const timestep +++122(t=3025, dt=25)\n",
+      "flux_top (\t3050.0\t)=\t-0.0040346727815483\n",
+      "flux_bot (\t3050.0\t)=\t0.0021388029068562194\n",
+      "mass_bot (\t3050.0\t)=\t2.7364105549159135\n",
+      "tlag=1770.5877894854254\n",
+      "+++ Const timestep +++123(t=3050, dt=25)\n",
+      "flux_top (\t3075.0\t)=\t-0.004020051847854768\n",
+      "flux_bot (\t3075.0\t)=\t0.0021529111014075026\n",
+      "mass_bot (\t3075.0\t)=\t2.79005698001921\n",
+      "tlag=1779.0537910761097\n",
+      "+++ Const timestep +++124(t=3075, dt=25)\n",
+      "flux_top (\t3100.0\t)=\t-0.004005663271857898\n",
+      "flux_bot (\t3100.0\t)=\t0.0021668162784223024\n",
+      "mass_bot (\t3100.0\t)=\t2.8440535722670828\n",
+      "tlag=1787.450523337452\n",
+      "+++ Const timestep +++125(t=3100, dt=25)\n",
+      "flux_top (\t3125.0\t)=\t-0.00399150291354155\n",
+      "flux_bot (\t3125.0\t)=\t0.002180520898995968\n",
+      "mass_bot (\t3125.0\t)=\t2.8983952869848113\n",
+      "tlag=1795.7784876909952\n",
+      "+++ Const timestep +++126(t=3125, dt=25)\n",
+      "flux_top (\t3150.0\t)=\t-0.00397756673127187\n",
+      "flux_bot (\t3150.0\t)=\t0.002194027422032807\n",
+      "mass_bot (\t3150.0\t)=\t2.953077140997671\n",
+      "tlag=1804.0381804974934\n",
+      "+++ Const timestep +++127(t=3150, dt=25)\n",
+      "flux_top (\t3175.0\t)=\t-0.003963850777952067\n",
+      "flux_bot (\t3175.0\t)=\t0.002207338302275277\n",
+      "mass_bot (\t3175.0\t)=\t3.008094212551522\n",
+      "tlag=1812.2300931620482\n",
+      "+++ Const timestep +++128(t=3175, dt=25)\n",
+      "flux_top (\t3200.0\t)=\t-0.0039503511979030225\n",
+      "flux_bot (\t3200.0\t)=\t0.002220455988489931\n",
+      "mass_bot (\t3200.0\t)=\t3.063441641186087\n",
+      "tlag=1820.3547122456382\n",
+      "+++ Const timestep +++129(t=3200, dt=25)\n",
+      "flux_top (\t3225.0\t)=\t-0.003937064223555061\n",
+      "flux_bot (\t3225.0\t)=\t0.0022333829217794096\n",
+      "mass_bot (\t3225.0\t)=\t3.119114627564454\n",
+      "tlag=1828.4125195694821\n",
+      "+++ Const timestep +++130(t=3225, dt=25)\n",
+      "flux_top (\t3250.0\t)=\t-0.003923986172721874\n",
+      "flux_bot (\t3250.0\t)=\t0.0022461215340219445\n",
+      "mass_bot (\t3250.0\t)=\t3.175108433261971\n",
+      "tlag=1836.4039923180085\n",
+      "+++ Const timestep +++131(t=3250, dt=25)\n",
+      "flux_top (\t3275.0\t)=\t-0.003911113445100562\n",
+      "flux_bot (\t3275.0\t)=\t0.0022586742464318\n",
+      "mass_bot (\t3275.0\t)=\t3.231418380517643\n",
+      "tlag=1844.3296031410634\n",
+      "+++ Const timestep +++132(t=3275, dt=25)\n",
+      "flux_top (\t3300.0\t)=\t-0.003898442520284095\n",
+      "flux_bot (\t3300.0\t)=\t0.0022710434682289117\n",
+      "mass_bot (\t3300.0\t)=\t3.2880398519509018\n",
+      "tlag=1852.1898202524933\n",
+      "+++ Const timestep +++133(t=3300, dt=25)\n",
+      "flux_top (\t3325.0\t)=\t-0.00388596995428224\n",
+      "flux_bot (\t3325.0\t)=\t0.0022832315954126737\n",
+      "mass_bot (\t3325.0\t)=\t3.3449682902464217\n",
+      "tlag=1859.9851075261383\n",
+      "+++ Const timestep +++134(t=3325, dt=25)\n",
+      "flux_top (\t3350.0\t)=\t-0.0038736923771592866\n",
+      "flux_bot (\t3350.0\t)=\t0.0022952410096370177\n",
+      "mass_bot (\t3350.0\t)=\t3.402199197809543\n",
+      "tlag=1867.715924591472\n",
+      "+++ Const timestep +++135(t=3350, dt=25)\n",
+      "flux_top (\t3375.0\t)=\t-0.0038616064917869034\n",
+      "flux_bot (\t3375.0\t)=\t0.0023070740771763904\n",
+      "mass_bot (\t3375.0\t)=\t3.4597281363947103\n",
+      "tlag=1875.3827269261142\n",
+      "+++ Const timestep +++136(t=3375, dt=25)\n",
+      "flux_top (\t3400.0\t)=\t-0.00384970907007086\n",
+      "flux_bot (\t3400.0\t)=\t0.002318733147974333\n",
+      "mass_bot (\t3400.0\t)=\t3.517550726709094\n",
+      "tlag=1882.9859659434899\n",
+      "+++ Const timestep +++137(t=3400, dt=25)\n",
+      "flux_top (\t3425.0\t)=\t-0.003837996951294282\n",
+      "flux_bot (\t3425.0\t)=\t0.0023302205547806607\n",
+      "mass_bot (\t3425.0\t)=\t3.5756626479935316\n",
+      "tlag=1890.52608908297\n",
+      "+++ Const timestep +++138(t=3425, dt=25)\n",
+      "flux_top (\t3450.0\t)=\t-0.003826467040222693\n",
+      "flux_bot (\t3450.0\t)=\t0.0023415386123618325\n",
+      "mass_bot (\t3450.0\t)=\t3.634059637582813\n",
+      "tlag=1898.003539895822\n",
+      "+++ Const timestep +++139(t=3450, dt=25)\n",
+      "flux_top (\t3475.0\t)=\t-0.003815116305072298\n",
+      "flux_bot (\t3475.0\t)=\t0.0023526896167803037\n",
+      "mass_bot (\t3475.0\t)=\t3.6927374904470898\n",
+      "tlag=1905.418758127277\n",
+      "+++ Const timestep +++140(t=3475, dt=25)\n",
+      "flux_top (\t3500.0\t)=\t-0.0038039417754467384\n",
+      "flux_bot (\t3500.0\t)=\t0.002363675844744135\n",
+      "mass_bot (\t3500.0\t)=\t3.751692058716145\n",
+      "tlag=1912.7721797985116\n",
+      "+++ Const timestep +++141(t=3500, dt=25)\n",
+      "flux_top (\t3525.0\t)=\t-0.0037929405405507724\n",
+      "flux_bot (\t3525.0\t)=\t0.002374499553021387\n",
+      "mass_bot (\t3525.0\t)=\t3.810919251188214\n",
+      "tlag=1920.0642372877762\n",
+      "+++ Const timestep +++142(t=3525, dt=25)\n",
+      "flux_top (\t3550.0\t)=\t-0.003782109747311109\n",
+      "flux_bot (\t3550.0\t)=\t0.0023851629779091136\n",
+      "mass_bot (\t3550.0\t)=\t3.870415032824845\n",
+      "tlag=1927.295359406536\n",
+      "+++ Const timestep +++143(t=3550, dt=25)\n",
+      "flux_top (\t3575.0\t)=\t-0.0037714465995350743\n",
+      "flux_bot (\t3575.0\t)=\t0.00239566833476023\n",
+      "mass_bot (\t3575.0\t)=\t3.930175424233212\n",
+      "tlag=1934.46597147532\n",
+      "+++ Const timestep +++144(t=3575, dt=25)\n",
+      "flux_top (\t3600.0\t)=\t-0.003760948354844343\n",
+      "flux_bot (\t3600.0\t)=\t0.0024060178175622495\n",
+      "mass_bot (\t3600.0\t)=\t3.990196501137243\n",
+      "tlag=1941.5764953976668\n",
+      "+++ Const timestep +++145(t=3600, dt=25)\n",
+      "flux_top (\t3625.0\t)=\t-0.003750612324323058\n",
+      "flux_bot (\t3625.0\t)=\t0.0024162135985617943\n",
+      "mass_bot (\t3625.0\t)=\t4.050474393838793\n",
+      "tlag=1948.6273497302714\n",
+      "+++ Const timestep +++146(t=3625, dt=25)\n",
+      "flux_top (\t3650.0\t)=\t-0.003740435871000156\n",
+      "flux_bot (\t3650.0\t)=\t0.0024262578279410285\n",
+      "mass_bot (\t3650.0\t)=\t4.111005286670078\n",
+      "tlag=1955.6189497557393\n",
+      "+++ Const timestep +++147(t=3650, dt=25)\n",
+      "flux_top (\t3675.0\t)=\t-0.0037304164073927422\n",
+      "flux_bot (\t3675.0\t)=\t0.0024361526335290587\n",
+      "mass_bot (\t3675.0\t)=\t4.171785417438454\n",
+      "tlag=1962.5517075483388\n",
+      "+++ Const timestep +++148(t=3675, dt=25)\n",
+      "flux_top (\t3700.0\t)=\t-0.003720551395736081\n",
+      "flux_bot (\t3700.0\t)=\t0.0024459001205577487\n",
+      "mass_bot (\t3700.0\t)=\t4.2328110768645395\n",
+      "tlag=1969.4260320411963\n",
+      "+++ Const timestep +++149(t=3700, dt=25)\n",
+      "flux_top (\t3725.0\t)=\t-0.003710838344964186\n",
+      "flux_bot (\t3725.0\t)=\t0.002455502371450594\n",
+      "mass_bot (\t3725.0\t)=\t4.294078608014644\n",
+      "tlag=1976.2423290888917\n",
+      "+++ Const timestep +++150(t=3725, dt=25)\n",
+      "flux_top (\t3750.0\t)=\t-0.003701274811171559\n",
+      "flux_bot (\t3750.0\t)=\t0.002464961445651902\n",
+      "mass_bot (\t3750.0\t)=\t4.355584405728425\n",
+      "tlag=1983.0010015322916\n",
+      "+++ Const timestep +++151(t=3750, dt=25)\n",
+      "flux_top (\t3775.0\t)=\t-0.0036918583952437484\n",
+      "flux_bot (\t3775.0\t)=\t0.002474279379482478\n",
+      "mass_bot (\t3775.0\t)=\t4.417324916042605\n",
+      "tlag=1989.7024492575547\n",
+      "+++ Const timestep +++152(t=3775, dt=25)\n",
+      "flux_top (\t3800.0\t)=\t-0.0036825867417521543\n",
+      "flux_bot (\t3800.0\t)=\t0.002483458186029655\n",
+      "mass_bot (\t3800.0\t)=\t4.479296635611507\n",
+      "tlag=1996.3470692564258\n",
+      "+++ Const timestep +++153(t=3800, dt=25)\n",
+      "flux_top (\t3825.0\t)=\t-0.00367345753860543\n",
+      "flux_bot (\t3825.0\t)=\t0.002492499855063495\n",
+      "mass_bot (\t3825.0\t)=\t4.541496111125172\n",
+      "tlag=2002.9352556835033\n",
+      "+++ Const timestep +++154(t=3825, dt=25)\n",
+      "flux_top (\t3850.0\t)=\t-0.0036644685155475046\n",
+      "flux_bot (\t3850.0\t)=\t0.0025014063529762275\n",
+      "mass_bot (\t3850.0\t)=\t4.603919938725668\n",
+      "tlag=2009.467399909725\n",
+      "+++ Const timestep +++155(t=3850, dt=25)\n",
+      "flux_top (\t3875.0\t)=\t-0.003655617442840878\n",
+      "flux_bot (\t3875.0\t)=\t0.002510179622754368\n",
+      "mass_bot (\t3875.0\t)=\t4.666564763422301\n",
+      "tlag=2015.9438905803179\n",
+      "+++ Const timestep +++156(t=3875, dt=25)\n",
+      "flux_top (\t3900.0\t)=\t-0.003646902130897349\n",
+      "flux_bot (\t3900.0\t)=\t0.00251882158396409\n",
+      "mass_bot (\t3900.0\t)=\t4.729427278506281\n",
+      "tlag=2022.365113664316\n",
+      "+++ Const timestep +++157(t=3900, dt=25)\n",
+      "flux_top (\t3925.0\t)=\t-0.0036383204285338807\n",
+      "flux_bot (\t3925.0\t)=\t0.002527334132763728\n",
+      "mass_bot (\t3925.0\t)=\t4.792504224965379\n",
+      "tlag=2028.7314525069908\n",
+      "+++ Const timestep +++158(t=3925, dt=25)\n",
+      "flux_top (\t3950.0\t)=\t-0.0036298702226187324\n",
+      "flux_bot (\t3950.0\t)=\t0.002535719141933438\n",
+      "mass_bot (\t3950.0\t)=\t4.855792390899094\n",
+      "tlag=2035.0432878790184\n",
+      "+++ Const timestep +++159(t=3950, dt=25)\n",
+      "flux_top (\t3975.0\t)=\t-0.0036215494378843984\n",
+      "flux_bot (\t3975.0\t)=\t0.002543978460919889\n",
+      "mass_bot (\t3975.0\t)=\t4.91928861093476\n",
+      "tlag=2041.3009980218262\n",
+      "+++ Const timestep +++160(t=3975, dt=25)\n",
+      "flux_top (\t4000.0\t)=\t-0.0036133560341165844\n",
+      "flux_bot (\t4000.0\t)=\t0.0025521139159087928\n",
+      "mass_bot (\t4000.0\t)=\t4.982989765645119\n",
+      "tlag=2047.5049586998136\n",
+      "+++ Const timestep +++161(t=4000, dt=25)\n",
+      "flux_top (\t4025.0\t)=\t-0.0036052880074166553\n",
+      "flux_bot (\t4025.0\t)=\t0.0025601273098981294\n",
+      "mass_bot (\t4025.0\t)=\t5.046892780967705\n",
+      "tlag=2053.6555432399464\n",
+      "+++ Const timestep +++162(t=4025, dt=25)\n",
+      "flux_top (\t4050.0\t)=\t-0.0035973433877324167\n",
+      "flux_bot (\t4050.0\t)=\t0.0025680204227959444\n",
+      "mass_bot (\t4050.0\t)=\t5.110994627626381\n",
+      "tlag=2059.7531225776775\n",
+      "+++ Const timestep +++163(t=4050, dt=25)\n",
+      "flux_top (\t4075.0\t)=\t-0.0035895202397565875\n",
+      "flux_bot (\t4075.0\t)=\t0.002575795011525168\n",
+      "mass_bot (\t4075.0\t)=\t5.1752923205553945\n",
+      "tlag=2065.7980652967317\n",
+      "+++ Const timestep +++164(t=4075, dt=25)\n",
+      "flux_top (\t4100.0\t)=\t-0.00358181666042476\n",
+      "flux_bot (\t4100.0\t)=\t0.0025834528101496747\n",
+      "mass_bot (\t4100.0\t)=\t5.23978291832633\n",
+      "tlag=2071.7907376745316\n",
+      "+++ Const timestep +++165(t=4100, dt=25)\n",
+      "flux_top (\t4125.0\t)=\t-0.00357423077950612\n",
+      "flux_bot (\t4125.0\t)=\t0.0025909955300015006\n",
+      "mass_bot (\t4125.0\t)=\t5.304463522578219\n",
+      "tlag=2077.731503718516\n",
+      "+++ Const timestep +++166(t=4125, dt=25)\n",
+      "flux_top (\t4150.0\t)=\t-0.0035667607582147097\n",
+      "flux_bot (\t4150.0\t)=\t0.0025984248598239783\n",
+      "mass_bot (\t4150.0\t)=\t5.369331277451038\n",
+      "tlag=2083.6207252055133\n",
+      "+++ Const timestep +++167(t=4150, dt=25)\n",
+      "flux_top (\t4175.0\t)=\t-0.0035594047891381437\n",
+      "flux_bot (\t4175.0\t)=\t0.0026057424659195224\n",
+      "mass_bot (\t4175.0\t)=\t5.434383369022831\n",
+      "tlag=2089.4587617160664\n",
+      "+++ Const timestep +++168(t=4175, dt=25)\n",
+      "flux_top (\t4200.0\t)=\t-0.0035521610950688216\n",
+      "flux_bot (\t4200.0\t)=\t0.0026129499923136976\n",
+      "mass_bot (\t4200.0\t)=\t5.499617024750747\n",
+      "tlag=2095.2459706735594\n",
+      "+++ Const timestep +++169(t=4200, dt=25)\n",
+      "flux_top (\t4225.0\t)=\t-0.0035450279282743433\n",
+      "flux_bot (\t4225.0\t)=\t0.002620049060920459\n",
+      "mass_bot (\t4225.0\t)=\t5.5650295129161735\n",
+      "tlag=2100.982707376822\n",
+      "+++ Const timestep +++170(t=4225, dt=25)\n",
+      "flux_top (\t4250.0\t)=\t-0.0035380035704574626\n",
+      "flux_bot (\t4250.0\t)=\t0.0026270412717185523\n",
+      "mass_bot (\t4250.0\t)=\t5.6306181420741614\n",
+      "tlag=2106.669325034724\n",
+      "+++ Const timestep +++171(t=4250, dt=25)\n",
+      "flux_top (\t4275.0\t)=\t-0.0035310863315147897\n",
+      "flux_bot (\t4275.0\t)=\t0.0026339282029298433\n",
+      "mass_bot (\t4275.0\t)=\t5.696380260507266\n",
+      "tlag=2112.3061747959146\n",
+      "+++ Const timestep +++172(t=4275, dt=25)\n",
+      "flux_top (\t4300.0\t)=\t-0.003524274549701532\n",
+      "flux_bot (\t4300.0\t)=\t0.002640711411212208\n",
+      "mass_bot (\t4300.0\t)=\t5.762313255684042\n",
+      "tlag=2117.8936057844826\n",
+      "+++ Const timestep +++173(t=4300, dt=25)\n",
+      "flux_top (\t4325.0\t)=\t-0.003517566590924664\n",
+      "flux_bot (\t4325.0\t)=\t0.002647392431843819\n",
+      "mass_bot (\t4325.0\t)=\t5.828414553722242\n",
+      "tlag=2123.43196512315\n",
+      "+++ Const timestep +++174(t=4325, dt=25)\n",
+      "flux_top (\t4350.0\t)=\t-0.0035109608467926794\n",
+      "flux_bot (\t4350.0\t)=\t0.002653972778928544\n",
+      "mass_bot (\t4350.0\t)=\t5.894681618856897\n",
+      "tlag=2128.9215979688065\n",
+      "+++ Const timestep +++175(t=4350, dt=25)\n",
+      "flux_top (\t4375.0\t)=\t-0.0035044557364769524\n",
+      "flux_bot (\t4375.0\t)=\t0.0026604539455928225\n",
+      "mass_bot (\t4375.0\t)=\t5.961111952913414\n",
+      "tlag=2134.3628475364894\n",
+      "+++ Const timestep +++176(t=4375, dt=25)\n",
+      "flux_top (\t4400.0\t)=\t-0.0034980497042966714\n",
+      "flux_bot (\t4400.0\t)=\t0.0026668374041914937\n",
+      "mass_bot (\t4400.0\t)=\t6.027703094785718\n",
+      "tlag=2139.7560551266006\n",
+      "+++ Const timestep +++177(t=4400, dt=25)\n",
+      "flux_top (\t4425.0\t)=\t-0.0034917412197734486\n",
+      "flux_bot (\t4425.0\t)=\t0.002673124606517718\n",
+      "mass_bot (\t4425.0\t)=\t6.094452619919583\n",
+      "tlag=2145.1015601517984\n",
+      "+++ Const timestep +++178(t=4425, dt=25)\n",
+      "flux_top (\t4450.0\t)=\t-0.0034855287780674662\n",
+      "flux_bot (\t4450.0\t)=\t0.002679316984012342\n",
+      "mass_bot (\t4450.0\t)=\t6.161358139801209\n",
+      "tlag=2150.3997001599914\n",
+      "+++ Const timestep +++179(t=4450, dt=25)\n",
+      "flux_top (\t4475.0\t)=\t-0.003479410898157669\n",
+      "flux_bot (\t4475.0\t)=\t0.0026854159479823554\n",
+      "mass_bot (\t4475.0\t)=\t6.228417301451143\n",
+      "tlag=2155.650810861995\n",
+      "+++ Const timestep +++180(t=4475, dt=25)\n",
+      "flux_top (\t4500.0\t)=\t-0.0034733861224642256\n",
+      "flux_bot (\t4500.0\t)=\t0.0026914228898116586\n",
+      "mass_bot (\t4500.0\t)=\t6.295627786923568\n",
+      "tlag=2160.8552261498653\n",
+      "flux_top (\t4525.0\t)=\t-0.0034674530176724227\n",
+      "flux_bot (\t4525.0\t)=\t0.0026973391811829154\n",
+      "mass_bot (\t4525.0\t)=\t6.362987312811\n",
+      "tlag=2166.0132781222874\n",
+      "+++ Const timestep +++181(t=4500, dt=25)\n",
+      "+++ Const timestep +++182(t=4525, dt=25)\n",
+      "flux_top (\t4550.0\t)=\t-0.003461610173329517\n",
+      "flux_bot (\t4550.0\t)=\t0.0027031661742946867\n",
+      "mass_bot (\t4550.0\t)=\t6.43049362975447\n",
+      "tlag=2171.1252971037484\n",
+      "+++ Const timestep +++183(t=4550, dt=25)\n",
+      "flux_top (\t4575.0\t)=\t-0.003455856201226943\n",
+      "flux_bot (\t4575.0\t)=\t0.002708905202084346\n",
+      "mass_bot (\t4575.0\t)=\t6.498144521959208\n",
+      "tlag=2176.1916116668604\n",
+      "+++ Const timestep +++184(t=4575, dt=25)\n",
+      "flux_top (\t4600.0\t)=\t-0.0034501897356454133\n",
+      "flux_bot (\t4600.0\t)=\t0.002714557578446086\n",
+      "mass_bot (\t4600.0\t)=\t6.565937806715839\n",
+      "tlag=2181.2125486487466\n",
+      "flux_top (\t4625.0\t)=\t-0.0034446094326966814+++ Const timestep +++185(t=4600, dt=25)\n",
+      "\n",
+      "flux_bot (\t4625.0\t)=\t0.0027201245984596846\n",
+      "mass_bot (\t4625.0\t)=\t6.633871333927161\n",
+      "tlag=2186.1884331755614\n",
+      "+++ Const timestep +++186(t=4625, dt=25)\n",
+      "flux_top (\t4650.0\t)=\t-0.0034391139696914297\n",
+      "flux_bot (\t4650.0\t)=\t0.00272560753860754\n",
+      "mass_bot (\t4650.0\t)=\t6.701942985640501\n",
+      "tlag=2191.1195886754876\n",
+      "+++ Const timestep +++187(t=4650, dt=25)\n",
+      "flux_top (\t4675.0\t)=\t-0.0034337020448653153\n",
+      "flux_bot (\t4675.0\t)=\t0.0027310076569974397\n",
+      "mass_bot (\t4675.0\t)=\t6.770150675585563\n",
+      "tlag=2196.006336895851\n",
+      "+++ Const timestep +++188(t=4675, dt=25)\n",
+      "flux_top (\t4700.0\t)=\t-0.0034283723775957466\n",
+      "flux_bot (\t4700.0\t)=\t0.0027363261935896906\n",
+      "mass_bot (\t4700.0\t)=\t6.838492348717902\n",
+      "tlag=2200.848997923481\n",
+      "+++ Const timestep +++189(t=4700, dt=25)\n",
+      "flux_top (\t4725.0\t)=\t-0.003423123706710885\n",
+      "flux_bot (\t4725.0\t)=\t0.0027415643704132254\n",
+      "mass_bot (\t4725.0\t)=\t6.906965980767938\n",
+      "tlag=2205.647890194576\n",
+      "+++ Const timestep +++190(t=4725, dt=25)\n",
+      "flux_top (\t4750.0\t)=\t-0.0034179547918561483\n",
+      "flux_bot (\t4750.0\t)=\t0.002746723391799483\n",
+      "mass_bot (\t4750.0\t)=\t6.975569577795597\n",
+      "tlag=2210.4033305204293\n",
+      "+++ Const timestep +++191(t=4750, dt=25)\n",
+      "flux_top (\t4775.0\t)=\t-0.0034128644116283562\n",
+      "flux_bot (\t4775.0\t)=\t0.0027518044445896486\n",
+      "mass_bot (\t4775.0\t)=\t7.044301175750461\n",
+      "tlag=2215.1156340886305\n",
+      "flux_top (\t4800.0\t)=\t-0.0034078513640181138\n",
+      "flux_bot (\t4800.0\t)=\t0.0027568086983670156\n",
+      "mass_bot (\t4800.0\t)=\t7.1131588400374195\n",
+      "tlag=2219.7851144873166\n",
+      "+++ Const timestep +++192(t=4775, dt=25)\n",
+      "+++ Const timestep +++193(t=4800, dt=25)\n",
+      "flux_top (\t4825.0\t)=\t-0.003402914466144534\n",
+      "flux_bot (\t4825.0\t)=\t0.0027617373056709744\n",
+      "mass_bot (\t4825.0\t)=\t7.1821406650878945\n",
+      "tlag=2224.4120837126593\n",
+      "+++ Const timestep +++194(t=4825, dt=25)\n",
+      "flux_top (\t4850.0\t)=\t-0.0033980525534890514\n",
+      "flux_bot (\t4850.0\t)=\t0.0027665914022213813\n",
+      "mass_bot (\t4850.0\t)=\t7.251244773936549\n",
+      "tlag=2228.9968521863034\n",
+      "+++ Const timestep +++195(t=4850, dt=25)\n",
+      "flux_top (\t4875.0\t)=\t-0.0033932644796267567\n",
+      "flux_bot (\t4875.0\t)=\t0.0027713721071267404\n",
+      "mass_bot (\t4875.0\t)=\t7.320469317803401\n",
+      "tlag=2233.5397287580404\n",
+      "+++ Const timestep +++196(t=4875, dt=25)\n",
+      "flux_top (\t4900.0\t)=\t-0.0033885491162762714\n",
+      "flux_bot (\t4900.0\t)=\t0.002776080523114475\n",
+      "mass_bot (\t4900.0\t)=\t7.389812475681416\n",
+      "tlag=2238.041020729899\n",
+      "+++ Const timestep +++197(t=4900, dt=25)\n",
+      "flux_top (\t4925.0\t)=\t-0.0033839053527894283\n",
+      "flux_bot (\t4925.0\t)=\t0.00278071773674293\n",
+      "mass_bot (\t4925.0\t)=\t7.459272453929634\n",
+      "tlag=2242.501033863753\n",
+      "+++ Const timestep +++198(t=4925, dt=25)\n",
+      "flux_top (\t4950.0\t)=\t-0.0033793320959881142\n",
+      "flux_bot (\t4950.0\t)=\t0.0027852848186042853\n",
+      "mass_bot (\t4950.0\t)=\t7.528847485871474\n",
+      "tlag=2246.9200723808913\n",
+      "+++ Const timestep +++199(t=4950, dt=25)\n",
+      "flux_top (\t4975.0\t)=\t-0.003374828269228641\n",
+      "flux_bot (\t4975.0\t)=\t0.0027897828235510583\n",
+      "mass_bot (\t4975.0\t)=\t7.598535831398415\n",
+      "tlag=2251.2984389851567\n",
+      "+++ Const timestep +++200(t=4975, dt=25)\n",
+      "flux_top (\t5000.0\t)=\t-0.003370392813181112\n",
+      "flux_bot (\t5000.0\t)=\t0.002794212790900305\n",
+      "mass_bot (\t5000.0\t)=\t7.668335776579057\n",
+      "tlag=2255.636434865688\n"
+     ]
+    }
+   ],
    "source": [
     "# Solve problem.\n",
     "try:\n",
@@ -520,9 +1586,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 52,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAG0CAYAAAActAwdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABA8UlEQVR4nO3dd3hUVeLG8XfSC8kklARCBkJvoTdB7FgQxS4CAiIWXBAUXQu6PxXR4LrroquiWBBFig3UtcIqoNIJgdA7hBpaMqmTZOb8/ohkiaBCSHJnku/neeYxM3Nn8nIImdd7zz3XZowxAgAA8EJ+VgcAAAD4PRQVAADgtSgqAADAa1FUAACA16KoAAAAr0VRAQAAXouiAgAAvBZFBQAAeC2KCgAA8FoUFQAA4LUsLSpFRUV68skn1ahRI4WGhqpx48YaP368PB6PlbEAAICXCLDym7/wwgt64403NG3aNLVp00YrV67UsGHDZLfbNWbMmD99vcfj0f79+xURESGbzVYJiQEAwLkyxigrK0txcXHy8/vjfSaWFpUlS5bouuuuU9++fSVJCQkJmjlzplauXHlGr9+/f78cDkdFRgQAABUkLS1N8fHxf7iNpUWlV69eeuONN7RlyxY1b95ca9as0c8//6xJkyaddnuXyyWXy1Vy/8SFn9PS0hQZGVkZkQEAwDlyOp1yOByKiIj4020tLSqPPvqoMjMz1bJlS/n7+8vtduu5557TgAEDTrt9UlKSnnnmmVMej4yMpKgAAOBjzmTahqWTaWfPnq3p06drxowZSk5O1rRp0/SPf/xD06ZNO+32jz/+uDIzM0tuaWlplZwYAABUJps5cfzEAg6HQ4899phGjhxZ8tiECRM0ffp0bdq06U9f73Q6ZbfblZmZyR4VAAB8xNl8flu6RyU3N/eU2b7+/v6cngwAACRZPEfl2muv1XPPPacGDRqoTZs2Wr16tV566SXdeeedVsYCAABewtJDP1lZWfrb3/6mOXPmKD09XXFxcRowYID+7//+T0FBQX/6eg79AADge87m89vSonKuKCoAAPgen5mjAgAA8EcoKgAAwGtRVAAAgNeiqAAAAK9FUQEAAF6LogIAALwWRQUAAJzCGKP/bjwkj8faVUwoKgAA4BQfr9qr4dNW6p4PVllaVigqAACglG3pWXrq8/WSpI4NouTnZ7MsC0UFAACUyC90a+SHq5VX6FavprV130VNLM1DUQEAACXG/2eDNh/KUu0aQXqpf3tL96ZIFBUAAPCr/6zdrxnL9kiSXrq1g2IiQixORFEBAACSdh/N0WOfpkqS/nJxE13YvI7FiYpRVAAAqOZcRW6NmrFa2a4idWkYrbGXN7c6UgmKCgAA1dzEbzYpdV+mosIC9cqAjgrw95564D1JAABApftu/UFN/WWXJOkfN7dXXFSotYF+g6ICAEA1tfd4rv768RpJ0l29Gql361iLE52KogIAQDVU6PZo9MzVcuYXqX28XY9c1dLqSKdFUQEAoBr65/dblLwnQxEhAXp1YCcFBXhnJfDOVAAAoML8uDldbyzcLkl64aZ2ctQMszjR76OoAABQjRzMzNdDHxXPSxl8XkNd3baexYn+GEUFAIBqosjt0ehZq3Usp0Ct60Xqib6trI70pygqAABUE6/8sE3Ldx5TeJC/Xh3YUSGB/lZH+lMUFQAAqoHF247o3z9slSQ9f2NbNa5Tw+JEZ4aiAgBAFXc4y6Uxs1NkjNS/i0PXdahvdaQzRlEBAKAK83iMxn6UosNZLjWLqaGn+7WxOtJZoagAAFCFTV64XT9tPaKQQD+9NqiTQoO8f17KySgqAABUUct3HtM/v98sSRrfL1HNYyMsTnT2KCoAAFRBx3IKNHrmanmMdH2HON3SJd7qSGVCUQEAoIoxxujhj9fooDNfjWuHa8INbWWz2ayOVSYUFQAAqph3ft6pHzalKyjAT/8e2FE1ggOsjlRmFBUAAKqQlLQMTfxmkyTpb9e0Vps4u8WJzg1FBQCAKiIzr1CjZiSryGN0ddu6ur17A6sjnTNLi0pCQoJsNtspt5EjR1oZCwAAn2OM0aOfrNXe43ly1AxV0o3tfHZeysksPWi1YsUKud3ukvvr1q3T5ZdfrltuucXCVAAA+J4Plu7Wt+sPKtDfplcHdJI9NNDqSOXC0qJSp06dUvcnTpyoJk2a6KKLLrIoEQAAvid1b6Ym/GejJOnRq1qqvSPK2kDlyGumARcUFGj69OkaO3bs7+6qcrlccrlcJfedTmdlxQMAwCtl5hXqLzNWqcDt0eWtYzW8VyOrI5Urr5lMO3fuXGVkZOiOO+743W2SkpJkt9tLbg6Ho/ICAgDgZYwxeuSTNUo7lqf46FD94+b2VWJeyslsxhhjdQhJuvLKKxUUFKQvv/zyd7c53R4Vh8OhzMxMRUZGVkZMAAC8xrs/79T4/2xQoL9Nn4zo6TOHfJxOp+x2+xl9fnvFoZ/du3dr/vz5+uyzz/5wu+DgYAUHB1dSKgAAvFdKWoaSvimelzLu6lY+U1LOllcc+pk6dapiYmLUt29fq6MAAOD1MnMLNfLDZBW6jfok1tUdPROsjlRhLC8qHo9HU6dO1dChQxUQ4BU7eAAA8FrGGD38yRrty8hTg5pheuHmqrFeyu+xvKjMnz9fe/bs0Z133ml1FAAAvN47P+/UvA2HFOTvp9cHdVJkSNVYL+X3WL4L44orrpCXzOcFAMCrJe85ftJ1fFopsb5vX8fnTFi+RwUAAPy54zkFGvVh8XV8rmlXT7ef19DqSJWCogIAgJfzeIwe+niN9mfmq1HtcCXd2LZKz0s5GUUFAAAvN+WnHfphU7qCAvz06sCOiqji81JORlEBAMCLrdh1TC9+t1mS9PS1bdQmrurPSzkZRQUAAC91NNul+2eslttjdF2HOA3oVv0uHUNRAQDAC3k8Rg9+tEYHnflqXCdcz99QfealnIyiAgCAF5q8cLsWbTms4IDi9VLCgy1fUcQSFBUAALzM0h1H9c/vi+eljL+ujVrWrb4X3qWoAADgRQ5nuTR65mp5jHRjx/q6tUv1m5dyMooKAABewu0xenB2itKzXGoaU0MTbkislvNSTkZRAQDAS7z24zb9vO2IQgP99fqgTgoLqp7zUk5GUQEAwAss3n5Ek+ZvkSQ9e32imsdGWJzIO1BUAACwWHpWvkbPTJHHSLd0jtfNneOtjuQ1KCoAAFjI7TEaMzNFR7JdahEbofHXJVodyatQVAAAsNDL/92qJTuOKizIX68N6qTQIH+rI3kVigoAABZZsDld//5hqyTp+RvaqmlMDYsTeR+KCgAAFtiXkacHZ6fIGGlQ9wa6vmN9qyN5JYoKAACVrKDIo5EfJut4bqHa1rfrb9e0tjqS16KoAABQyZ7/eqNS0jIUGRKg1wd1Ukgg81J+D0UFAIBK9J+1+/Xe4l2SpJdu7SBHzTBrA3k5igoAAJVk++FsPfrJWknSfRc3Ue/WsRYn8n4UFQAAKkFuQZHum75KOQVudW9UUw9d3tzqSD6BogIAQAUzxujJOeu05VC26kQE698DOyrAn4/gM8EoAQBQwWatSNNnq/fJzyb9e0BHxUSEWB3JZ1BUAACoQOv2ZeqpL9ZLkv56ZUud17iWxYl8C0UFAIAKkplbqPs+XKWCIo96t4rRvRc2tjqSz6GoAABQAYwxeujjNUo7lqf46FD985YO8vOzWR3L51BUAACoAFMW7dD8jYcU5O+nyYM6yx4WaHUkn0RRAQCgnC3bcVR//26zJOn/rm2ttvF2ixP5LooKAADlKD0rX6NmrpbbY3R9hzgN6t7A6kg+jaICAEA5KXJ7NGZmig5nudQspoaeu6GtbDbmpZwLigoAAOXkX/O3aMmOowoL8tfk2zspPDjA6kg+j6ICAEA5+GHTIb3243ZJ0sSb2qlpTITFiaoGy4vKvn37dPvtt6tWrVoKCwtThw4dtGrVKqtjAQBwxtKO5erB2WskSUN6NFS/9nEWJ6o6LN0ndfz4cZ1//vm65JJL9M033ygmJkbbt29XVFSUlbEAADhjriK3Rs5IVmZeodrH2/VE31ZWR6pSLC0qL7zwghwOh6ZOnVryWEJCgnWBAAA4S8/+Z4PW7s2UPTRQrw3qpOAAf6sjVSmWHvr54osv1KVLF91yyy2KiYlRx44d9dZbb/3u9i6XS06ns9QNAACrfLpqr6Yv3SNJmtS/g+KjwyxOVPVYWlR27NihyZMnq1mzZvruu+80YsQIjR49Wu+///5pt09KSpLdbi+5ORyOSk4MAECxDfudGjcnVZI0+rJmuqRljMWJqiabMcZY9c2DgoLUpUsXLV68uOSx0aNHa8WKFVqyZMkp27tcLrlcrpL7TqdTDodDmZmZioyMrJTMAABk5hWq36s/a/fRXF3UvI7evaOr/LmOzxlzOp2y2+1n9Plt6R6VevXqqXXr1qUea9Wqlfbs2XPa7YODgxUZGVnqBgBAZfJ4jB76KEW7j+aqflSoJvXvQEmpQJYWlfPPP1+bN28u9diWLVvUsGFDixIBAPDHXl+wTfM3pisowE9v3N5Z0eFBVkeq0iwtKg8++KCWLl2q559/Xtu2bdOMGTM0ZcoUjRw50spYAACc1k9bD+uf87ZIkp69rg0XG6wElhaVrl27as6cOZo5c6YSExP17LPPatKkSRo0aJCVsQAAOMW+jDyNnrlaxkj9uzjUvysXG6wMlk6mPVdnMxkHAICychW5desbS7Rmb6YS60fqkxE9FRLIeill5TOTaQEA8AXPfLlBa35d1G3yoM6UlEpEUQEA4A98smqvZizbI5tNmnRbBzlqsqhbZaKoAADwO9bvz9QTvy7qNuayZrqkBYu6VTaKCgAAp5GZW6j7pifLVeTRxS3qaPSlzayOVC1RVAAA+A2Px+jBj1K051iu4qOLF3XzY1E3S1BUAAD4jdd+3KYfNv1vUbeoMBZ1swpFBQCAkyzaclgvzS9e1G3CdYlKrM+iblaiqAAA8Ku9x3M1Zlbxom4Dujl0a1eH1ZGqPYoKAACS8gvd+suHyTqeW6i29e166to2VkeCKCoAAEgqXtRt7d5MRYUFavLtnVjUzUtQVAAA1d7HK9M0c3nxom4v39ZR8dEs6uYtKCoAgGpt3b5MPTl3nSTpwd7NdVHzOhYnwskoKgCAautYToHu/WCVXEUeXdKijkZd0tTqSPgNigoAoFoqcns0euZq7cvIU8NaYZrUvyOLunkhigoAoFp68fvN+nnbEYUG+mvK4C6yhwVaHQmnQVEBAFQ7X609oDcX7pAkvXhLO7WoG2FxIvweigoAoFrZfDBLf/1kjSTpngsb65p2cRYnwh+hqAAAqo3MvELd+8FK5Ra4dX7TWnrkyhZWR8KfoKgAAKoFj8do7OwU7Tqaq/pRofr3gE4K8Odj0NvxNwQAqBZe/u9W/XdTuoID/PTm4M6qGc4VkX0BRQUAUOXN23BIL/93qyTp+RvackVkH0JRAQBUadsPZ2vs7BRJ0tAeDXVT53hrA+GsUFQAAFVWtqtI936wSlmuInVNiNaT17S2OhLOEkUFAFAlGWP014/XaFt6tmIjg/XaoE4KZPKsz+FvDABQJU1euF3frDuoQH+bXh/UWTERIVZHQhlQVAAAVc6iLYf1j+82S5Ke7tdGnRtGW5wIZUVRAQBUKWnHcnX/zNXyGKl/F4cGdmtgdSScA4oKAKDKyCtw654PVikzr1Dt4+165ro2stm4IrIvo6gAAKoEY4we/2ytNh5wqlZ4kCbf3lkhgf5Wx8I5oqgAAKqEqb/s0tyU/fL3s+m1QZ0UFxVqdSSUA4oKAMDnLdl+VM99vVGS9MTVrXRe41oWJ0J5oagAAHza/ow8jZqRLLfH6PoOcRp2foLVkVCOLC0qTz/9tGw2W6lb3bp1rYwEAPAhxZNnV+poToFa1YtU0o3tmDxbxQRYHaBNmzaaP39+yX1/fyY+AQD+nDFGj322Vuv2OVUzPEhTBndWaBCfIVWN5UUlICCAvSgAgLP21k879PmJybMDO8lRM8zqSKgAls9R2bp1q+Li4tSoUSPddttt2rFjh9WRAABebuGWw5r4zSZJ0v9d01o9mjB5tqqydI9K9+7d9f7776t58+Y6dOiQJkyYoJ49e2r9+vWqVevUHzqXyyWXy1Vy3+l0VmZcAIAX2HUkR/fPSJbHSLd2ideQHg2tjoQKZOkelT59+uimm25S27Zt1bt3b3311VeSpGnTpp12+6SkJNnt9pKbw+GozLgAAItlu4p09/sr5cwvUscGUXr2+kQmz1Zxlh/6OVl4eLjatm2rrVu3nvb5xx9/XJmZmSW3tLS0Sk4IALCKx2M0dnaKtqZnKyYiWG/c3lnBAUyereosn0x7MpfLpY0bN+qCCy447fPBwcEKDg6u5FQAAG/wyg9b9f2GQwry99ObgzsrNjLE6kioBJbuUXn44Ye1cOFC7dy5U8uWLdPNN98sp9OpoUOHWhkLAOBlvlt/UJPmF+9tn3BDojo2iLY4ESqLpXtU9u7dqwEDBujIkSOqU6eOzjvvPC1dulQNGzIxCgBQbMuhLI2dnSJJuqNngm7twvzE6sTSojJr1iwrvz0AwMtl5Bbo7vdXKqfArR6Na+mJvq2sjoRK5lWTaQEAOKHI7dH9M1dr99Fc1Y8K1WuDOinQn4+t6oa/cQCAV/r7d5v109YjCg3015QhnVUzPMjqSLAARQUA4HXmrt6nKYuKVyp/8ZZ2ahNntzgRrEJRAQB4ldS9mXr007WSpL9c3ETXtIuzOBGsRFEBAHiNI9ku3fvBSrmKPLqkRR09dEULqyPBYhQVAIBXKCjy6L7pq7Q/M1+Na4fr5QEd5e/H8vjVHUUFAGA5Y4yenJuqFbuOKyI4QFOGdFFkSKDVseAFKCoAAMu9+8sufbRyr/xs0r8HdlTTmBpWR4KXoKgAACy1cMthPffVBknSuKtb6eIWMRYngjehqAAALLP9cLZGzUiWx0i3dI7X8F6NrI4EL0NRAQBYIjO3UHdNW6ms/CJ1aRitCTckymZj8ixKo6gAACpdkdujkTOStfNIjupHheqNwZ0VHOBvdSx4IYoKAKDSTfhqo37e9r/l8WvXCLY6ErwURQUAUKlmLt+j9xbvkiT9q397lsfHH6KoAAAqzbIdR/W3ueskSWMvb66rEutZnAjejqICAKgUacdydd+HySryGF3Trp7uv7Sp1ZHgAygqAIAKl+0q0l3TVupYToHa1rfrxZvbc4YPzghFBQBQoTweowdnp2jzoSzViQjWlCGdFRrEGT44MxQVAECF+ue8zZq34ZCCAvw0ZXBn1bOHWh0JPoSiAgCoMJ+n7NNrP26XJL1wU1t1bBBtcSL4GooKAKBCpKRl6K+frJUkjbioiW7oGG9xIvgiigoAoNwdzMzXPe+vVEGRR71bxeivV7awOhJ8FEUFAFCucguKNHzaCqVnudQ8tob+1b+D/P04wwdlQ1EBAJQbj8fogVkpWr/fqVrhQXpnaFdFhARaHQs+jKICACg3L3y3Sd+fOMNnSGc5aoZZHQk+jqICACgXH61I05sLd0iSXry5nTo3rGlxIlQFFBUAwDlbsv2oxs1JlSSNvqyZrutQ3+JEqCooKgCAc7LzSI7u+3BVyTV8HuzdzOpIqEIoKgCAMsvILdDw91YoI7dQHRxR+sctXMMH5YuiAgAok0K3R/dNT9aOIzmqHxWqKUM6KySQa/igfJWpqCQnJys1NbXk/ueff67rr79e48aNU0FBQbmFAwB4J2OM/jZ3nZbsOKrwIH+9PbSLYiJCrI6FKqhMReXee+/Vli1bJEk7duzQbbfdprCwMH388cd65JFHyjUgAMD7vP3TTs1akSY/m/TvgR3Vql6k1ZFQRZWpqGzZskUdOnSQJH388ce68MILNWPGDL333nv69NNPyzMfAMDLzNtwSM9/s1GS9GTf1rq0ZazFiVCVlamoGGPk8XgkSfPnz9fVV18tSXI4HDpy5EiZgiQlJclms+mBBx4o0+sBABVv/f5MjZm1WsZIg7o30LDzE6yOhCquTEWlS5cumjBhgj744AMtXLhQffv2lSTt3LlTsbFn36xXrFihKVOmqF27dmWJAwCoBOnOfN01baVyC9y6oFltPd2vDWf4oMKVqahMmjRJycnJGjVqlJ544gk1bdpUkvTJJ5+oZ8+eZ/Ve2dnZGjRokN566y1FR0eXJQ4AoILlFbh11/srdSAzX03qhOvVgZ0U6M+Jo6h4AWV5Ubt27Uqd9XPCiy++KH//szs1beTIkerbt6969+6tCRMmlCUOAKACeTxGD32corV7MxUdFqh37+gqeygXGkTlKFNRSUtLk81mU3x8vCRp+fLlmjFjhlq3bq177rnnjN9n1qxZSk5O1ooVK85oe5fLJZfLVXLf6XSeXXAAwFn757zN+jr1oAL9bXpzcBc1rBVudSRUI2Xabzdw4ED9+OOPkqSDBw/q8ssv1/LlyzVu3DiNHz/+jN4jLS1NY8aM0fTp0xUScmbn3iclJclut5fcHA5HWeIDAM7Q7BV79NqP2yVJE29sp26NuNAgKpfNGGPO9kXR0dFaunSpWrRooVdeeUWzZ8/WL7/8ou+//14jRozQjh07/vQ95s6dqxtuuKHUoSK32y2bzSY/Pz+5XK5TDiOdbo+Kw+FQZmamIiM5hx8AytNPWw/rjqkr5PYYjb6smcZe3tzqSKginE6n7Hb7GX1+l+nQT2FhoYKDgyUVn57cr18/SVLLli114MCBM3qPyy677JR5LsOGDVPLli316KOPnnauS3BwcMn3BQBUnM0Hs/SX6clye4xu6FifCw3CMmUqKm3atNEbb7yhvn37at68eXr22WclSfv371etWrXO6D0iIiKUmJhY6rHw8HDVqlXrlMcBAJXnkDNfw6YuV5arSN0b1dTEm9pyGjIsU6Y5Ki+88ILefPNNXXzxxRowYIDat28vSfriiy/UrVu3cg0IAKg8Oa4iDZ+2Qvsz89W4TrimDO6i4AAuNAjrlGmOilQ8n8TpdJZa+2TXrl0KCwtTTExMuQX8I2dzjAsA8MfcHqN73l+p/25KV63wIM35y/lqUCvM6liogip8jook+fv7n7JAW0JCQlnfDgBgIWOMxn+5Xv/dlK7gAD+9NbQLJQVeocxF5ZNPPtFHH32kPXv2qKCgoNRzycnJ5xwMAFB53v1ll6Yt2S2bTZrUv4M6NWClcHiHMs1ReeWVVzRs2DDFxMRo9erV6tatm2rVqqUdO3aoT58+5Z0RAFCBvl13UBO+2iBJGtenlfq0rWdxIuB/ylRUXn/9dU2ZMkWvvvqqgoKC9Mgjj2jevHkaPXq0MjMzyzsjAKCCpKRl6IHZxVdDvv28BrrrgkZWRwJKKVNR2bNnT8nFB0NDQ5WVlSVJGjx4sGbOnFl+6QAAFSbtWK7umrZC+YUeXdKijp6+lqshw/uUqajUrVtXR48elSQ1bNhQS5culSTt3LlTZTyJCABQiTJzCzXsvRU6kl2g1vUi9e+BnRTA1ZDhhcr0U3nppZfqyy+/lCQNHz5cDz74oC6//HL1799fN9xwQ7kGBACUL1eRWyOmr9K29GzVjQzRu3d0VY3gMp9bAVSoMq2j4vF45PF4FBBQ/IP98ccf66efflLTpk113333KTCwci7/zToqAHB2PB6jB2an6Is1+1UjOEAf3dtDreP4/YnKdTaf32Ve8C0/P19r165Venq6PB7P/97QZtO1115blrc8axQVADg7Sd9s1JsLdyjAz6apw7rqgmZ1rI6EaqjCF3z79ttvNXjw4JJ5Kiez2Wxyu91leVsAQAWatniX3lxYfHX7F25qR0mBTyjTHJVRo0bp1ltv1YEDB0oOA524UVIAwPt8u+6Anv5yvSTpr1e20E2d4y1OBJyZMhWV9PR0jR07VrGxseWdBwBQzlbuOqYxs1JkjDSoewP95eImVkcCzliZisrNN9+sBQsWlHMUAEB525aerbveXylXkUe9W8XomX6slQLfUqbJtLm5ubrllltUp04dtW3b9pSzfEaPHl1uAf8Ik2kB4PelZ+XrxtcXa+/xPHVwRGnm3ecpNMjf6lhAxU+mnTFjhr777juFhoZqwYIFpdq5zWartKICADi9bFeRhk1dob3H85RQK0zvDO1CSYFPKlNRefLJJzV+/Hg99thj8vNjJUMA8CaFbo/+8mGy1u93qlZ4kKbd2U21agRbHQsokzK1jIKCAvXv35+SAgBexhijxz9L1aIthxUa6K937+iqhrXCrY4FlFmZmsbQoUM1e/bs8s4CADhH/5q3RZ+s2is/m/TaoI5q74iyOhJwTsp06Mftduvvf/+7vvvuO7Vr1+6UybQvvfRSuYQDAJy5Gcv26JUftkmSnruhrS5tyRIS8H1lKiqpqanq2LGjJGndunWlnuO0NwCofPM2HNKTc1MlSaMvbaoB3RpYnAgoH2UqKj/++GN55wAAlNGKXcc0akayPEa6uXO8Hry8udWRgHLDbFgA8GGbDjo1/L0VchV5dFnLGE28sS17tlGlUFQAwEftPZ6roe8ulzO/SJ0bRuvVgZ0U4M+vdVQt/EQDgA86llOgIe8u1yGnS81ja7CgG6osigoA+JgcV5GGTV2uHYdzFGcP0bQ7uykqLMjqWECFoKgAgA8pKPJoxPRVWrM3U1FhgXp/eDfVs4daHQuoMBQVAPARHo/RXz9Zo5+2HlFooL+m3tFVTWMirI4FVCiKCgD4AGOMJny1UZ+n7FeAn02Tb++kjg2irY4FVDiKCgD4gDcW7tC7v+yUJL14Sztd3CLG4kRA5aCoAICX+2hlml74dpMk6cm+rXRDx3iLEwGVh6ICAF5s/oZDevyz4qXx772ose66oLHFiYDKRVEBAC+1YtcxjZyRLLfH6KZO8XrsqpZWRwIqHUUFALzQ+v2ZuvPXpfEvbRmjiTexND6qJ0uLyuTJk9WuXTtFRkYqMjJSPXr00DfffGNlJACw3I7D2RryznJl5RepW0JNvTawkwJZGh/VlKU/+fHx8Zo4caJWrlyplStX6tJLL9V1112n9evXWxkLACyzLyNPt7+9TEdzCtQmLlJv38HS+KjebMYYY3WIk9WsWVMvvviihg8f/qfbOp1O2e12ZWZmKjIyshLSAUDFOZLt0q1vLNGOIzlqXCdcH9/bQ7VqBFsdCyh3Z/P5HVBJmf6U2+3Wxx9/rJycHPXo0eO027hcLrlcrpL7TqezsuIBQIXKzCvUkHeWa8eRHNWPCtX04d0pKYC8YDJtamqqatSooeDgYI0YMUJz5sxR69atT7ttUlKS7HZ7yc3hcFRyWgAof3kFbt01bYU2HHCqdo0gTb+ru+KiuH4PIHnBoZ+CggLt2bNHGRkZ+vTTT/X2229r4cKFpy0rp9uj4nA4OPQDwGcVFHl09/srtXDLYUWEBGj2PT3UOo7fZ6jazubQj+VF5bd69+6tJk2a6M033/zTbZmjAsCXuT1Go2eu1lepBxQa6K/pd3VT54Y1rY4FVLiz+fy2/NDPbxljSu01AYCqyBijJ+ak6qvUAwr0t+nNwZ0pKcBpWDqZdty4cerTp48cDoeysrI0a9YsLViwQN9++62VsQCgQhljlPTNJs1akSY/m/TKbR11YfM6VscCvJKlReXQoUMaPHiwDhw4ILvdrnbt2unbb7/V5ZdfbmUsAKhQry/YrimLdkiSJt7YTn3a1rM4EeC9LC0q77zzjpXfHgAq3QdLdunF7zZLKr4S8q1dOXsR+CNeN0cFAKqqj1ak6W+fF6+8PfqyZlwJGTgDFBUAqASfp+zTo5+tlSQN79VID/ZuZnEiwDdQVACggn2TekBjP1ojY6Tbz2ugJ/u24krIwBmiqABABfrvxkO6f+ZquT1Gt3SO1/h+iZQU4CxQVACggvy09bDum56sIo9Rv/ZxmnhTO/n5UVKAs0FRAYAKsHTHUd39/koVuD26sk2s/nlre/lTUoCzRlEBgHK2avdxDX9vhfILPbqkRR39e0AnBfrz6xYoC/7lAEA5St2bqTveXa6cArfOb1pLk2/vrKAAftUCZcW/HgAoJxsPODX43WXKchWpW0JNvTWki0IC/a2OBfg0igoAlINt6dm6/e1lysgtVAdHlN4d1lVhQZYu/g1UCRQVADhHu47kaOBbS3U0p0Bt4iI17c5uqhFMSQHKA0UFAM5B2rFcDXp7mdKzXGoRG6EPhneXPTTQ6lhAlUFRAYAySjuWq9umLNW+jDw1rh2uD+7qpprhQVbHAqoUigoAlMFvS8rMe85TTESI1bGAKoeiAgBn6XQlJTaSkgJUBIoKAJyFk0tKI0oKUOEoKgBwhvYez9WAt04qKXdTUoCKRlEBgDOw93jxnpS9x/OUUCtMM+8+T3XtlBSgolFUAOBP7MvI04C3Tiop91BSgMpCUQGAP7AvI0+3TVmitGN5avhrSalnD7U6FlBtsHQiAPyOE3NSTpSUWZQUoNJRVADgNHYfzdHAt5ZpX0aeGtQsnpNCSQEqH0UFAH5j++FsDXxrqQ45XWpcO1wzmDgLWIaiAgAn2XIoSwPfWqYj2S41i6mhD+/uzoqzgIUoKgDwqw37nbr9nWU6llOgVvUiNX14N9WqEWx1LKBao6gAgKTUvZm6/Z1lyswrVNv6dn0wvJuiwrjAIGA1igqAai95z3ENfXe5svKL1LFBlN4b1k320ECrYwEQRQVANbdi1zHd8e5y5RS41TUhWlOHdVONYH41At6Cf40Aqq3F249o+HsrlVfoVo/GtfTOHV0UFsSvRcCb8C8SQLX046Z0jZi+Sq4ijy5sXkdTBndWSKC/1bEA/AZFBUC18+Wa/XpwdoqKPEa9W8Xo1YGdKCmAl6KoAKhWZi3fo8fnpMoYqV/7OP3z1vYK9OeyZ4C3oqgAqDbeWrRDz329UZI0sHsDPXtdovz9bBanAvBHLP3fiKSkJHXt2lURERGKiYnR9ddfr82bN1sZCUAVZIzRS99vLikp917UWM9dT0kBfIGlRWXhwoUaOXKkli5dqnnz5qmoqEhXXHGFcnJyrIwFoArxeIye+XKDXvlhmyTpr1e20GNXtZTNRkkBfIHNGGOsDnHC4cOHFRMTo4ULF+rCCy/80+2dTqfsdrsyMzMVGRlZCQkB+JIit0ePfZaqT1btlSSNv66NhvRIsDYUgLP6/PaqOSqZmZmSpJo1a572eZfLJZfLVXLf6XRWSi4AvsdV5NYDs1L0zbqD8vez6cWb2+nGTvFWxwJwlrxmqrsxRmPHjlWvXr2UmJh42m2SkpJkt9tLbg6Ho5JTAvAFOa4i3TVtpb5Zd1BB/n56bWAnSgrgo7zm0M/IkSP11Vdf6eeff1Z8/Ol/oZxuj4rD4eDQD4ASx3IKNOy9FVqTlqHQQH+9NaSLejWrbXUsACfxuUM/999/v7744gstWrTod0uKJAUHBys4mEuuAzi9vcdzNeSd5dpxJEdRYYF6946u6tQg2upYAM6BpUXFGKP7779fc+bM0YIFC9SoUSMr4wDwYZsOOjX03eU65HQpzh6i94d3V9OYGlbHAnCOLC0qI0eO1IwZM/T5558rIiJCBw8elCTZ7XaFhoZaGQ2AD1m+85iGT1uhrPwiNY+toWl3dlM9O79DgKrA0jkqv7eOwdSpU3XHHXf86es5PRnAd+sP6v6Zq1VQ5FGXhtF6Z2hX2cMCrY4F4A/4zBwVL5nHC8BHzVq+R+PmpMpjpN6tYvXqwI5cXBCoYrxiMi0AnA1jjF79YZv+OW+LJKl/F4eeuyFRAVxcEKhyKCoAfIrbYzT+y/WatmS3JGnUJU310BXNWRIfqKIoKgB8Rl6BW6Nnrda8DYdks0lPXdNad5zP2YJAVUZRAeATDme5dNe0FVqzN1NBAX566db2uqZdnNWxAFQwigoAr7ctPVvD3luutGN5ig4L1FtDuqhLwumvCQagaqGoAPBqy3Yc1T0frFJmXqEa1grTe8O6qVHtcKtjAagkFBUAXuvzlH3668drVeD2qGODKL09pItq1eAyGkB1QlEB4HWMMXp9wXa9+N1mSVKfxLr6V/8OrJECVEMUFQBepdDt0d/mrtOsFWmSpLt6NdK4q1vJz4/Tj4HqiKICwGs48ws1asZqLdpyWH426alr22hozwSrYwGwEEUFgFfYfTRHw6et1Lb0bIUG+uvfAzqqd+tYq2MBsBhFBYDllmw/qvs+XKWM3ELFRgbr7SFd1TbebnUsAF6AogLAUjOX79Hf5q5TkceofbxdU4Z0UWxkiNWxAHgJigoASxS5PXru642a+ssuSdK17eP04s3tOLMHQCkUFQCV7uRJs5L00OXNNerSplxYEMApKCoAKtWuIzkaPm2Fth/OUWigv166tb36tK1ndSwAXoqiAqDSLN5+RPdNT1ZmXqHq2UP01pAuSqzPpFkAv4+iAqDCGWP03uJdeu6rjSryGHVwRGnK4M6KYdIsgD9BUQFQofIK3Hr8s7Wam7JfknR9hzhNvIlJswDODEUFQIVJO5arez9YpQ0HnPL3s2nc1a105/kJTJoFcMYoKgAqxKIth3X/zNXKzCtUrfAgvTqwk3o0qWV1LAA+hqICoFyduPLxP77fLGOk9o4ovXF7J9Wzh1odDYAPoqgAKDfZriI9/NEafbv+oCTptq4OPd2vDfNRAJQZRQVAudh+OFv3vL9S2w/nKNDfpmf6JWpg9wZWxwLg4ygqAM7Z16kH9Mgna5XtKlJsZLAm395ZnRpEWx0LQBVAUQFQZq4it57/aqOmLdktSeqWUFOvDuqomAjWRwFQPigqAMok7ViuRs5I1tq9mZKkERc10cNXNFeAv5/FyQBUJRQVAGdt3oZDeuijFDnzi2QPDdS/+rfXpS1jrY4FoAqiqAA4Y4Vuj/7+7Sa99dNOSVLHBlF6dWAn1Y/i1GMAFYOiAuCMpB3L1ZhZq5W8J0OSNLxXIz16VUsFBXCoB0DFoagA+FNfrNmvJz5LVZarSBEhAXrx5va6KrGu1bEAVAMUFQC/K8dVpKe+WK9PVu2VJHVuGK1J/TvIUTPM4mQAqguKCoDTSt2bqdGzVmvnkRz52aRRlzbT6EubclYPgEpl6W+cRYsW6dprr1VcXJxsNpvmzp1rZRwAkjweoymLtuvGyb9o55Ec1bOHaObd52ns5Zx6DKDyWfpbJycnR+3bt9err75qZQwAv9qfkach7y7X819vUqHb6Ko2dfXNmAvUvTFXPQZgDUsP/fTp00d9+vSxMgIAFV/xeG7KPv3f5+uVlV+kkEA//d81bTSgm0M2m83qeACqMZ+ao+JyueRyuUruO51OC9MAVcOxnAI9OTdVX6cWX/G4gyNKL93aXo3r1LA4GQD4WFFJSkrSM888Y3UMoMr4YdMhPfJJqo5kuxTgZ9OYy5rpvoubMBcFgNfwqaLy+OOPa+zYsSX3nU6nHA6HhYkA35TtKtJzX23QzOVpkqSmMTX0r1s7qG283eJkAFCaTxWV4OBgBQcHWx0D8GkLtxzWuM9StS8jTzabNPz8Rnr4yhYKCfS3OhoAnMKnigqAssvMLdSzX20oWbzNUTNUL9zUTj2b1LY4GQD8PkuLSnZ2trZt21Zyf+fOnUpJSVHNmjXVoEEDC5MBVct36w/qybnrdDjLJZtNuqNngv56ZQuFBfH/KgC8m6W/pVauXKlLLrmk5P6J+SdDhw7Ve++9Z1EqoOo4mu3SU1+s13/WHpAkNa4Trr/f1E5dEmpanAwAzoylReXiiy+WMcbKCECVZIzRJ6v26vmvN+p4bqH8/Wy658LGGnNZM+aiAPAp7PcFqpht6VkaN2edlu88JklqWTdCL97cnjN6APgkigpQReQXuvXqD9v05qLtKnQbhQT66YHezTW8VyMFsi4KAB9FUQGqgIVbDutvc9dpz7FcSdJlLWP0dL82ctQMszgZAJwbigrgw/Zl5On5rzbqq9TiybL17CF66to2urJNLNfoAVAlUFQAH5RX4Nabi7brjYXblV/okb+fTcN6JuiBy5urRjD/rAFUHfxGA3yIMUZfpR5Q0tebtC8jT5LUvVFNPXVtG7WOi7Q4HQCUP4oK4CM27HfqmS/Xa9mvZ/PUjwrVuKtb6eq2dTnMA6DKoqgAXu6QM1+T5m/R7BVp8hgpJNBPIy5qonsvbKLQINZEAVC1UVQAL5WVX6gpi3bo7Z92Kq/QLUnq266exl3dSvWjQi1OBwCVg6ICeJmCIo9mLNutV37YpmM5BZKkTg2i9PjVrdSVpe8BVDMUFcBLnJgo++J3m7X7aPF6KI1rh+uRq1pyujGAaouiAljMGKN5Gw5p0vyt2nDAKUmqXSNYD/Rupv5dHawqC6Bao6gAFjHG6IdN6Zo0f6tS92VKksKD/HX3hY119wWNFc56KABAUQEqmzFGC7Yc1qR5W7Rmb3FBCQvy19CeCbr7gsaqGR5kcUIA8B4UFaCSeDxG/92UrtcXbNPqPRmSpNBAfw3p2VD3XNBYtWoEWxsQALwQRQWoYAVFHn2xZr/eXLhdW9OzJRWvhTL4vIa696Imqk1BAYDfRVEBKkiOq0gzl+/ROz/v1IHMfElSRHCABp3XUHf2SlBMRIjFCQHA+1FUgHJ2MDNfHy7brfeX7FZmXqEkqU5EsO48v5EGnddAkSGBFicEAN9BUQHKgTFGq3Yf19TFu/TduoMq8hhJUkKtMN17URPd0LG+QgJZ7h4AzhZFBTgH+YVufbFmv6Yt3qX1+50lj3dNiNaw8xvpyjZ15e/HQm0AUFYUFaAMtqVn6+OVafpoZZqO5xYf3gkO8NN1HeI0tGeC2sTZLU4IAFUDRQU4QzmuIn2VekAfrUjTyt3HSx6vHxWq289rqNu6OhTNGigAUK4oKsAfMMZodVqGPlqRpi/X7FdOQfFVjP1s0iUtYnRrV4cuaxmjAJa5B4AKQVEBTmProSx9sWa/vlyzX7t+vUCgVDw59pYuDt3cOV6xkZxeDAAVjaIC/CrtWG5JOdl0MKvk8ZBAP12dWE+3dnWoe6OaXMUYACoRRQXVljFG2w9na96GdH23/qBS0jJKngv0t+nCZnXUr0OcereK5QKBAGARfvuiWnF7itc7mbfhoOZvTNfOIzklz9lsUo/GtdSvfZyuSqyrqDAmxgKA1SgqqPIOOfP1y7Yj+nnrEf24Ob3kdGJJCvL3U48mtdS7dayubB2rGOadAIBXoaigysl2FWn5zqP6aesR/bLtiLYcyi71vD00UJe2jNHlrWN1QbPaimBJewDwWhQV+LzDWS6t2n1cq3Yf06rdx7V2b2bJEvZS8SGdxDi7ejWrrQub1VHXhGhOJwYAH0FRgU9xFbm19VC21u7N1Mrdx5S8+3ip04dPcNQMVa+mddSraW31bFKLhdgAwEdRVOC1cguKtPGAU+v2ObV+f6bW7XNqa3qWCt2m1HY2m9Q8JkKdGkarS8NodU2oqQa1wixKDQAoTxQVWO5otkvbD+doW3q2tqVna/vh4v/uy8g77fb20EAl1o9UpwbR6tQwWp0aRMseyjwTAKiKLC8qr7/+ul588UUdOHBAbdq00aRJk3TBBRdYHQvlxBijLFeRDmTka+/xXKUdy9Xe43naezxPab/ed+YX/e7rYyOD1SbOrsS4SLWOsyuxfqTqR4Wy6BoAVBOWFpXZs2frgQce0Ouvv67zzz9fb775pvr06aMNGzaoQYMGVkbD7zDGKNtVpIzcQmXmFd9OfH08t0DpznylZ7l+veXrcJZL+YWeP3xPm634wn5NY2qoSZ0apf5bk7klAFCt2Ywx5s83qxjdu3dXp06dNHny5JLHWrVqpeuvv15JSUl/+nqn0ym73a7MzExFRkaWW67cgiIdyynQb0fGGMnIlHwtSabkOXPS1yWvOGXb/339m+dOevx03/e3rzn5+5bOUfxfjzEqLPKo0GNU5Pao0G1U6PaoyFP8ddGv9wt/fS6voEi5BW7lFrqVX+Au+brk8QJ3STFxe87+RyYyJECOmmFyRIcpPjpU8dGhctQMU3x0mBrUDFNokP9ZvycAwDedzee3ZXtUCgoKtGrVKj322GOlHr/iiiu0ePHi077G5XLJ5XKV3Hc6nRWSbf7GdI2eubpC3ruqCArwU1RooKLCAmUPDZQ9NEhRYYGKiQguvkWG/Pp1iOpEBFNEAABlYllROXLkiNxut2JjY0s9Hhsbq4MHD572NUlJSXrmmWcqPJu/zaaQwOJ1NmwqngtxYkrEiZkRNput5Gud5rk/e43tNy8+ddsT9//3/G/fX7/d9qTX+NmkQH8/Bfj7KdDfVvy1X/F/A/1tpR4P9PdTaKC/woL8FRr0638D/RUaFPC/xwL9FRUWJPuv5SQkkOIBAKh4lk+m/e2HrjHmdydKPv744xo7dmzJfafTKYfDUe6Z+rarp77t6pX7+wIAgLNjWVGpXbu2/P39T9l7kp6efspelhOCg4MVHBxcGfEAAIAXsGwd8aCgIHXu3Fnz5s0r9fi8efPUs2dPi1IBAABvYumhn7Fjx2rw4MHq0qWLevTooSlTpmjPnj0aMWKElbEAAICXsLSo9O/fX0ePHtX48eN14MABJSYm6uuvv1bDhg2tjAUAALyEpeuonKuKWkcFAABUnLP5/OZa9wAAwGtRVAAAgNeiqAAAAK9FUQEAAF6LogIAALwWRQUAAHgtigoAAPBaFBUAAOC1KCoAAMBrWbqE/rk6saiu0+m0OAkAADhTJz63z2RxfJ8uKllZWZIkh8NhcRIAAHC2srKyZLfb/3Abn77Wj8fj0f79+xURESGbzVau7+10OuVwOJSWlsZ1hCoQ41w5GOfKwThXDsa58lTUWBtjlJWVpbi4OPn5/fEsFJ/eo+Ln56f4+PgK/R6RkZH8Q6gEjHPlYJwrB+NcORjnylMRY/1ne1JOYDItAADwWhQVAADgtSgqvyM4OFhPPfWUgoODrY5SpTHOlYNxrhyMc+VgnCuPN4y1T0+mBQAAVRt7VAAAgNeiqAAAAK9FUQEAAF6LogIAALwWReU0Xn/9dTVq1EghISHq3LmzfvrpJ6sjebVFixbp2muvVVxcnGw2m+bOnVvqeWOMnn76acXFxSk0NFQXX3yx1q9fX2obl8ul+++/X7Vr11Z4eLj69eunvXv3ltrm+PHjGjx4sOx2u+x2uwYPHqyMjIwK/tN5h6SkJHXt2lURERGKiYnR9ddfr82bN5fahnEuH5MnT1a7du1KFrjq0aOHvvnmm5LnGefyl5SUJJvNpgceeKDkMca5fDz99NOy2WylbnXr1i153ifG2aCUWbNmmcDAQPPWW2+ZDRs2mDFjxpjw8HCze/duq6N5ra+//to88cQT5tNPPzWSzJw5c0o9P3HiRBMREWE+/fRTk5qaavr372/q1atnnE5nyTYjRoww9evXN/PmzTPJycnmkksuMe3btzdFRUUl21x11VUmMTHRLF682CxevNgkJiaaa665prL+mJa68sorzdSpU826detMSkqK6du3r2nQoIHJzs4u2YZxLh9ffPGF+eqrr8zmzZvN5s2bzbhx40xgYKBZt26dMYZxLm/Lly83CQkJpl27dmbMmDEljzPO5eOpp54ybdq0MQcOHCi5paenlzzvC+NMUfmNbt26mREjRpR6rGXLluaxxx6zKJFv+W1R8Xg8pm7dumbixIklj+Xn5xu73W7eeOMNY4wxGRkZJjAw0MyaNatkm3379hk/Pz/z7bffGmOM2bBhg5Fkli5dWrLNkiVLjCSzadOmCv5TeZ/09HQjySxcuNAYwzhXtOjoaPP2228zzuUsKyvLNGvWzMybN89cdNFFJUWFcS4/Tz31lGnfvv1pn/OVcebQz0kKCgq0atUqXXHFFaUev+KKK7R48WKLUvm2nTt36uDBg6XGNDg4WBdddFHJmK5atUqFhYWltomLi1NiYmLJNkuWLJHdblf37t1LtjnvvPNkt9ur5d9NZmamJKlmzZqSGOeK4na7NWvWLOXk5KhHjx6MczkbOXKk+vbtq969e5d6nHEuX1u3blVcXJwaNWqk2267TTt27JDkO+Ps0xclLG9HjhyR2+1WbGxsqcdjY2N18OBBi1L5thPjdrox3b17d8k2QUFBio6OPmWbE68/ePCgYmJiTnn/mJiYavd3Y4zR2LFj1atXLyUmJkpinMtbamqqevToofz8fNWoUUNz5sxR69atS37pMs7nbtasWUpOTtaKFStOeY6f5/LTvXt3vf/++2revLkOHTqkCRMmqGfPnlq/fr3PjDNF5TRsNlup+8aYUx7D2SnLmP52m9NtXx3/bkaNGqW1a9fq559/PuU5xrl8tGjRQikpKcrIyNCnn36qoUOHauHChSXPM87nJi0tTWPGjNH333+vkJCQ392OcT53ffr0Kfm6bdu26tGjh5o0aaJp06bpvPPOk+T948yhn5PUrl1b/v7+pzTA9PT0UxonzsyJ2eV/NKZ169ZVQUGBjh8//ofbHDp06JT3P3z4cLX6u7n//vv1xRdf6Mcff1R8fHzJ44xz+QoKClLTpk3VpUsXJSUlqX379nr55ZcZ53KyatUqpaenq3PnzgoICFBAQIAWLlyoV155RQEBASVjwDiXv/DwcLVt21Zbt271mZ9nispJgoKC1LlzZ82bN6/U4/PmzVPPnj0tSuXbGjVqpLp165Ya04KCAi1cuLBkTDt37qzAwMBS2xw4cEDr1q0r2aZHjx7KzMzU8uXLS7ZZtmyZMjMzq8XfjTFGo0aN0meffaYffvhBjRo1KvU841yxjDFyuVyMczm57LLLlJqaqpSUlJJbly5dNGjQIKWkpKhx48aMcwVxuVzauHGj6tWr5zs/z+c8HbeKOXF68jvvvGM2bNhgHnjgARMeHm527dpldTSvlZWVZVavXm1Wr15tJJmXXnrJrF69uuSU7okTJxq73W4+++wzk5qaagYMGHDa09/i4+PN/PnzTXJysrn00ktPe/pbu3btzJIlS8ySJUtM27Ztq81phvfdd5+x2+1mwYIFpU4zzM3NLdmGcS4fjz/+uFm0aJHZuXOnWbt2rRk3bpzx8/Mz33//vTGGca4oJ5/1YwzjXF4eeughs2DBArNjxw6zdOlSc80115iIiIiSzzRfGGeKymm89tprpmHDhiYoKMh06tSp5BRQnN6PP/5oJJ1yGzp0qDGm+BS4p556ytStW9cEBwebCy+80KSmppZ6j7y8PDNq1ChTs2ZNExoaaq655hqzZ8+eUtscPXrUDBo0yERERJiIiAgzaNAgc/z48Ur6U1rrdOMryUydOrVkG8a5fNx5550l//7r1KljLrvsspKSYgzjXFF+W1QY5/JxYl2UwMBAExcXZ2688Uazfv36kud9YZxtxhhz7vtlAAAAyh9zVAAAgNeiqAAAAK9FUQEAAF6LogIAALwWRQUAAHgtigoAAPBaFBUAAOC1KCoAKt2CBQtks9mUkZFhdRQAXo4F3wBUuIsvvlgdOnTQpEmTJBVfT+TYsWOKjY2tNlexBVA2AVYHAFD9BAUFlVy5FQD+CId+AFSoO+64QwsXLtTLL78sm80mm82m9957r9Shn/fee09RUVH6z3/+oxYtWigsLEw333yzcnJyNG3aNCUkJCg6Olr333+/3G53yXsXFBTokUceUf369RUeHq7u3btrwYIF1vxBAVQI9qgAqFAvv/yytmzZosTERI0fP16StH79+lO2y83N1SuvvKJZs2YpKytLN954o2688UZFRUXp66+/1o4dO3TTTTepV69e6t+/vyRp2LBh2rVrl2bNmqW4uDjNmTNHV111lVJTU9WsWbNK/XMCqBgUFQAVym63KygoSGFhYSWHezZt2nTKdoWFhZo8ebKaNGkiSbr55pv1wQcf6NChQ6pRo4Zat26tSy65RD/++KP69++v7du3a+bMmdq7d6/i4uIkSQ8//LC+/fZbTZ06Vc8//3zl/SEBVBiKCgCvEBYWVlJSJCk2NlYJCQmqUaNGqcfS09MlScnJyTLGqHnz5qXex+VyqVatWpUTGkCFo6gA8AqBgYGl7ttsttM+5vF4JEkej0f+/v5atWqV/P39S213crkB4NsoKgAqXFBQUKlJsOWhY8eOcrvdSk9P1wUXXFCu7w3Ae3DWD4AKl5CQoGXLlmnXrl06cuRIyV6Rc9G8eXMNGjRIQ4YM0WeffaadO3dqxYoVeuGFF/T111+XQ2oA3oCiAqDCPfzww/L391fr1q1Vp04d7dmzp1zed+rUqRoyZIgeeughtWjRQv369dOyZcvkcDjK5f0BWI+VaQEAgNdijwoAAPBaFBUAAOC1KCoAAMBrUVQAAIDXoqgAAACvRVEBAABei6ICAAC8FkUFAAB4LYoKAADwWhQVAADgtSgqAADAa1FUAACA1/p/7z1I6n0j4D8AAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
     "import matplotlib.pyplot as pyplot\n",
     "pyplot.xlabel(\"time\")\n",
@@ -569,7 +1646,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,
diff --git a/content/skin/SkinDiffusionWithLagtime.py b/content/skin/SkinDiffusionWithLagtime.py
index 18c86f4..6ebbba9 100644
--- a/content/skin/SkinDiffusionWithLagtime.py
+++ b/content/skin/SkinDiffusionWithLagtime.py
@@ -1,31 +1,38 @@
 #!/usr/bin/env python
 # coding: utf-8
 
-# [<img src="../header.svg">](../index.ipynd)
+# [<img src="../../header.svg">](../index.ipynd)
 # 
 # ---
 # 
 # # Drug Transport across a Virtual Skin Membrane:  Extended version
-# [Previous Version](SkinDiffusion.ipybnd)
+# We know the basic steps from the [previous version](SkinDiffusion.ipybnd).
 
 # ## 1. Setup
 # 
 # ### Load existing modules 
 
-# In[1]:
+# In[27]:
 
 
-import modsimtools as util
+import ug4py.pyugcore  as ugcore
+import ug4py.pylimex as limex
+import ug4py.pyconvectiondiffusion as cd
+# import pysuperlu as slu
+
+
+# In[28]:
 
-import ug4py as ug4
-import pylimex as limex
-import pyconvectiondiffusion as cd
-import pysuperlu as slu
+
+# Importing
+import sys
+sys.path.append('..')
+import modsimtools as util
 
 
 # ### Create Domain
 
-# In[2]:
+# In[29]:
 
 
 requiredSubsets = {"LIP", "COR", "BOTTOM_SC", "TOP_SC"}
@@ -33,7 +40,7 @@
 numRefs = 2
 
 
-# In[3]:
+# In[30]:
 
 
 dom = util.CreateDomain(gridName, numRefs, requiredSubsets)
@@ -41,17 +48,16 @@
 
 # ### Create Approximation space
 
-# In[4]:
+# In[31]:
 
 
 approxSpaceDesc = dict(fct = "u", type = "Lagrange", order = 1)
 approxSpace = util.CreateApproximationSpace(dom, approxSpaceDesc)
 
 
-# ### Create discretization
-# 
+# ### Create spatial discretization
 
-# In[5]:
+# In[32]:
 
 
 # Define model parameter
@@ -59,7 +65,7 @@
 D={ "COR":  1.0, "LIP":1.0 }
 
 
-# In[6]:
+# In[33]:
 
 
 # Create element discretizations for lipids and corneocytes.
@@ -68,20 +74,20 @@
 elemDisc["LIP"] = util.CreateDiffusionElemDisc("u", "LIP", K["LIP"], K["LIP"]*D["LIP"], 0.0)
 
 
-# In[7]:
+# In[34]:
 
 
 # Set Dirichlet boundary conditions.
-dirichletBnd = ug4.DirichletBoundary2dCPU1()
+dirichletBnd = ugcore.DirichletBoundary2dCPU1()
 dirichletBnd.add(1.0, "u", "TOP_SC")
 dirichletBnd.add(0.0, "u", "BOTTOM_SC")
 
 
-# In[8]:
+# In[35]:
 
 
-# Create the global discretization object.
-domainDisc = ug4.DomainDiscretization2dCPU1(approxSpace)
+# Create the global spatial discretization object.
+domainDisc = ugcore.DomainDiscretization2dCPU1(approxSpace)
 domainDisc.add(elemDisc["LIP"])
 domainDisc.add(elemDisc["COR"])
 domainDisc.add(dirichletBnd)
@@ -90,47 +96,47 @@
 # ## 2. Steady state problem
 # Flux is computed from steady state. Since configuration of a multigrid solver is somewhat tricky, we use an LU decomposition here:
 
-# In[9]:
+# In[36]:
 
 
-A = ug4.AssembledLinearOperatorCPU1(domainDisc)
-u = ug4.GridFunction2dCPU1(approxSpace)
-b = ug4.GridFunction2dCPU1(approxSpace)
+A = ugcore.AssembledLinearOperatorCPU1(domainDisc)
+u = ugcore.GridFunction2dCPU1(approxSpace)
+b = ugcore.GridFunction2dCPU1(approxSpace)
 
-ug4.Interpolate(0.0, u, "u")
+ugcore.Interpolate(0.0, u, "u")
 
 domainDisc.assemble_linear(A, b)
 domainDisc.adjust_solution(u)
 
-#lsolver = ug4.LUCPU1()
-lsolver =slu.SuperLUCPU1()
+lsolver = ugcore.LUCPU1()
+#lsolver =slu.SuperLUCPU1()
 
 lsolver.init(A, u)
 lsolver.apply(u, b)
 
-ug4.WriteGridFunctionToVTK(u, "SkinSteadyState.vtk")
+ugcore.WriteGridFunctionToVTK(u, "vtk/SkinSteadyState.vtk")
 
 
 # Compute $J_\infty=J(t=\infty)$ for
 # $$ J(t)=\frac{1}{|\Gamma|}\int_\Gamma (-KD \nabla u(t,x)) \cdot \vec n dA$$
 
-# In[10]:
+# In[37]:
 
 
-areaG=ug4.Integral(1.0, u, "BOTTOM_SC")
+areaG=ugcore.Integral(1.0, u, "BOTTOM_SC")
 print("Surface area [um^2]:")
 print(areaG)
 
 surfaceFlux = {}
-surfaceFlux["BOT"] = K["LIP"]*D["LIP"]*ug4.IntegrateNormalGradientOnManifold(u, "u", "BOTTOM_SC", "LIP")
-surfaceFlux["TOP"] = K["LIP"]*D["LIP"]*ug4.IntegrateNormalGradientOnManifold(u, "u", "TOP_SC", "LIP")
+surfaceFlux["BOT"] = K["LIP"]*D["LIP"]*ugcore.IntegrateNormalGradientOnManifold(u, "u", "BOTTOM_SC", "LIP")
+surfaceFlux["TOP"] = K["LIP"]*D["LIP"]*ugcore.IntegrateNormalGradientOnManifold(u, "u", "TOP_SC", "LIP")
 
 print("Surface fluxes [kg/s]:")
 print(surfaceFlux["TOP"])
 print(surfaceFlux["BOT"])
 
 
-# In[11]:
+# In[38]:
 
 
 print("Normalized Fluxes [kg / (mu^2 * s)]:")
@@ -138,7 +144,7 @@
 print(surfaceFlux["BOT"]/areaG)
 
 
-# In[12]:
+# In[39]:
 
 
 Jref = 1.0/17.6
@@ -150,19 +156,47 @@
 
 # ## 3. Transient problem
 
-# After each time-step, we execute a a callback function `MyPostProcess`. In this function, print the solution and compute
+# ### Define time integrator
+
+# In[40]:
+
+
+# Create time discretization.
+timeDisc=ugcore.ThetaTimeStepCPU1(domainDisc, 1.0) 
+
+
+# In[41]:
+
+
+# Settings for time stepping.
+startTime = 0.0
+endTime = 5000.0
+dt=25.0
+dtMin=2.5
+
+
+# In[42]:
+
+
+# Create time integrator.
+timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
+timeInt.set_linear_solver(lsolver)
+timeInt.set_time_step(dt)
+
+
+# ###  Creating observers
+
+# After each successful step, we would like to perform certain operations, e.g. for book-keeping. To that end, we would like to execute the following callback `MyPostProcess`. In this function, print the solution and compute
 # $$
 # m(t_k):= \int_0^{t_k} J(s) \, ds \approx \sum_{i=1}^k(t_{i}- t_{i-1}) \frac{J(t_{i-1}) +J(t_i)}{2} 
 # $$
 # using the trapezoid rule. Moreover, we also compute the lag time $\tau$ from $m(t_k) = J_\infty(t_k - \tau)$.
-# 
 
-# In[13]:
+# In[43]:
 
 
 # auxiliary variables for book-keeping
 import numpy as np
-import os
 
 tPoints = np.array([0.0])
 mPoints = np.array([0.0])
@@ -180,8 +214,8 @@ def MyPostProcess(u, step, time, dt):
     
     # 1) Compute fluxes.
     gradFlux={}
-    gradFlux["BOT"] = ug4.IntegrateNormalGradientOnManifold(u, "u", "BOTTOM_SC", "LIP")
-    gradFlux["TOP"] = ug4.IntegrateNormalGradientOnManifold(u, "u", "TOP_SC", "LIP")
+    gradFlux["BOT"] = ugcore.IntegrateNormalGradientOnManifold(u, "u", "BOTTOM_SC", "LIP")
+    gradFlux["TOP"] = ugcore.IntegrateNormalGradientOnManifold(u, "u", "TOP_SC", "LIP")
   
     jTOP = K["LIP"]*D["LIP"]*gradFlux["TOP"]
     jBOT = K["LIP"]*D["LIP"]*gradFlux["BOT"]
@@ -200,87 +234,82 @@ def MyPostProcess(u, step, time, dt):
     tPoints = np.append(tPoints, time)
     jPoints = np.append(jPoints, jBOT)
     mPoints = np.append(mPoints, mass)
-    
-    os.system('df -k /')
-    return True
-    
-
 
-# In[14]:
 
+# In[44]:
 
-pyobserver= ug4.PythonCallbackObserver2dCPU1(MyPostProcess)
 
+# We use this function to define a callback
+pyobserver= ugcore.PythonCallbackObserver2dCPU1(MyPostProcess)
 
-# Additionally, we define a callback for vtk files.
 
-# In[15]:
+# A second observer is used for writing vtk files:
 
+# In[45]:
 
-# Callback for file I/O.
-vtk=ug4.VTKOutput2d()
-vtkobserver=limex.VTKOutputObserver2dCPU1("vtk/SkinData.vtu", vtk)
 
+# Create callback observer for file I/O.
+def MyVTKCallback(usol, step, time, dt) :
+    ugcore.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_"+str(int(step)).zfill(5)+".vtu")
+    
+vtkobserver = ugcore.PythonCallbackObserver2dCPU1(MyVTKCallback) 
 
-# ### Define time integrator
-
-# In[16]:
 
+# In[46]:
 
-# Create time discretization.
-timeDisc=ug4.ThetaTimeStepCPU1(domainDisc, 1.0) 
 
+# Alternative:
+# vtk=ug4.VTKOutput2d()
+# vtkobserver=limex.VTKOutputObserver2dCPU1("../../SkinData.vtu", vtk)
 
-# In[17]:
 
+# The third  observer (for convenience, non-essential) indicates the progress of time stepping:
 
-# Settings for time stepping.
-startTime = 0.0
-endTime = 2500.0
-dt=25.0
-dtMin=2.5
+# In[47]:
 
-ug4.Interpolate(0.0, u, "u")
 
+# Create progress widget.
+import ipywidgets as widgets
+from IPython.display import display
 
-# In[18]:
+wProgress=widgets.FloatProgress(value=startTime, min=startTime, max=endTime, 
+                                description='Progress:', style={'bar_color': 'blue'})
 
+# Create callback and attach observer.
+def MyProgressCallback(u, step, time, dt):
+    wProgress.value=time
+    
+progressobserver = ugcore.PythonCallbackObserver2dCPU1(MyProgressCallback)
 
-# Create time integrator.
-timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
-timeInt.set_linear_solver(lsolver)
-timeInt.set_time_step(dt)
 
+# All observers must be attached to the time integrator:
 
-# In[19]:
+# In[48]:
 
 
 # Attach observers.
-timeInt.attach_observer(vtkobserver)
 timeInt.attach_observer(pyobserver)
+timeInt.attach_observer(vtkobserver)
+timeInt.attach_observer(progressobserver) 
 
 
-# In[20]:
+# ### Solve transient problem
 
+# In[49]:
 
-# Create progress widget.
-import ipywidgets as widgets
-from IPython.display import display
 
-wProgress=widgets.FloatProgress(value=startTime, min=startTime, max=endTime, 
-                                description='Progress:', style={'bar_color': 'blue'})
-display(wProgress)
+# Set initial values.
+ugcore.Interpolate(0.0, u, "u")
 
-# Create callback and attach observer.
-def MyProgress(u, step, time, dt):
-    wProgress.value=time
-    
-timeInt.attach_observer(ug4.PythonCallbackObserver2dCPU1(MyProgress)) 
 
+# In[50]:
 
-# ### Solve transient problem
 
-# In[ ]:
+# Show progress bar.
+display(wProgress)
+
+
+# In[51]:
 
 
 # Solve problem.
@@ -295,7 +324,7 @@ def MyProgress(u, step, time, dt):
 # ## 4. Visualization
 # 
 
-# In[ ]:
+# In[52]:
 
 
 import matplotlib.pyplot as pyplot

From a3382cf3356244e84ef9b8f5701fd5ab0e7c0d98 Mon Sep 17 00:00:00 2001
From: Arne Naegel <naegel@gcsc.uni-frankfurt.de>
Date: Fri, 5 Jul 2024 00:41:10 +0200
Subject: [PATCH 4/4] Some geometry files

---
 content/tutorial-fem/grids/sectorTest.ugx     | 33 +++++++++++++++++++
 content/tutorial-fem/grids/sector_ref0.ugx    | 33 +++++++++++++++++++
 .../tutorial-fem/grids/unit_square_quad.ugx   | 18 ++++++++++
 .../tutorial-fem/grids/unit_square_tri.ugx    | 29 ++++++++++++++++
 4 files changed, 113 insertions(+)
 create mode 100644 content/tutorial-fem/grids/sectorTest.ugx
 create mode 100644 content/tutorial-fem/grids/sector_ref0.ugx
 create mode 100644 content/tutorial-fem/grids/unit_square_quad.ugx
 create mode 100644 content/tutorial-fem/grids/unit_square_tri.ugx

diff --git a/content/tutorial-fem/grids/sectorTest.ugx b/content/tutorial-fem/grids/sectorTest.ugx
new file mode 100644
index 0000000..40eebbe
--- /dev/null
+++ b/content/tutorial-fem/grids/sectorTest.ugx
@@ -0,0 +1,33 @@
+<?xml version="1.0" encoding="utf-8"?>
+<grid name="defGrid">
+	<vertices coords="3">0 0 0 0 1 0 0.500000013513858255 0.866025395982208845 0 0.866025419388897655 0.499999972972283213 0 1 6.12323399573676604e-17 0 1.22464679914735321e-16 -1 0 -0.499999891889129133 -0.866025466202270056 0 -0.866025466202269723 -0.499999891889129799 0 -1 -1.83697019872102969e-16 0 -0.866025466202269945 0.499999891889129411 0 -0.499999891889129799 0.866025466202269723 0 0.433012709694448827 0.249999986486141607 0 0.5 3.06161699786838302e-17 0 0.250000006756929127 0.433012697991104423 0 0 0.5 0 -0.2499999459445649 0.433012733101134861 0 -0.433012733101134972 0.249999945944564705 0 -0.5 -9.18485099360514844e-17 0 -0.433012733101134861 -0.2499999459445649 0 -0.249999945944564567 -0.433012733101135028 0 6.12323399573676604e-17 -0.5 0 0.96592583032779944 0.258819030029770414 0 0.707106797737575898 0.707106764635518692 0 0.25881905263889593 0.965925824269702527 0 -0.258818984811518715 0.9659258424439916 0 -0.707106781186547684 0.707106781186547351 0 -0.965925842443991711 0.258818984811518438 0 -0.9659258424439916 -0.258818984811518826 0 -0.707106781186547351 -0.707106781186547906 0 -0.258818984811518327 -0.965925842443991822 0</vertices>
+	<edges>0 11 11 3 0 12 12 4 0 13 13 2 1 14 14 0 0 15 15 10 0 16 16 9 0 17 17 8 0 18 18 7 0 19 19 6 0 20 20 5 11 13 12 11 13 14 14 15 15 16 16 17 17 18 18 19 19 20 4 21 21 3 21 12 21 11 3 22 22 2 2 23 23 1 1 24 24 10 10 25 25 9 9 26 26 8 8 27 27 7 7 28 28 6 6 29 29 5 22 11 22 13 23 13 23 14 24 14 24 15 25 15 25 16 26 16 26 17 27 17 27 18 28 18 28 19 29 19 29 20</edges>
+	<triangles>13 0 11 11 0 12 14 0 13 15 0 14 16 0 15 17 0 16 18 0 17 19 0 18 20 0 19 12 4 21 12 21 11 11 21 3 11 3 22 11 22 13 13 22 2 13 2 23 13 23 14 14 23 1 14 1 24 14 24 15 15 24 10 15 10 25 15 25 16 16 25 9 16 9 26 16 26 17 17 26 8 17 8 27 17 27 18 18 27 7 18 7 28 18 28 19 19 28 6 19 6 29 19 29 20 20 29 5</triangles>
+	<subset_handler name="defSH">
+		<subset name="Inner" color="0.588235 0.588235 1 1" state="262144">
+			<vertices>11 13 14 15 16 17 18 19</vertices>
+			<edges>0 1 4 5 6 7 8 9 10 11 12 13 14 15 16 17 20 21 22 23 24 25 26 27 28 31 32 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64</edges>
+			<faces>0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35</faces>
+		</subset>
+		<subset name="Circle" color="1 0 0 1" state="262144">
+			<vertices>1 2 3 4 5 6 7 8 9 10 21 22 23 24 25 26 27 28 29</vertices>
+			<edges>29 30 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48</edges>
+		</subset>
+		<subset name="Cut" color="0 1 0 1" state="393216">
+			<vertices>0 12 20</vertices>
+			<edges>2 3 18 19</edges>
+		</subset>
+	</subset_handler>
+	<subset_handler name="markSH">
+		<subset name="crease" color="1 1 1 1" state="0"/>
+		<subset name="fixed" color="1 1 1 1" state="0"/>
+	</subset_handler>
+	<selector name="defSel">
+		<vertices>16 1</vertices>
+	</selector>
+	<projection_handler name="defPH" subset_handler="0">
+		<default type="default">0 0</default>
+		<projector type="sphere" subset="1">0 0 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 -1.00000000000000000e+00</projector>
+	</projection_handler>
+</grid>
+
diff --git a/content/tutorial-fem/grids/sector_ref0.ugx b/content/tutorial-fem/grids/sector_ref0.ugx
new file mode 100644
index 0000000..51ca0c7
--- /dev/null
+++ b/content/tutorial-fem/grids/sector_ref0.ugx
@@ -0,0 +1,33 @@
+<?xml version="1.0" encoding="utf-8"?>
+<grid name="defGrid">
+	<vertices coords="3">0 0 0 0 1 0 0.25881905263889593 0.965925824269702638 0 0.500000013513858255 0.866025395982208845 0 0.707106781186547462 0.707106781186547573 0 0.866025419388897655 0.499999972972283213 0 0.965925818211605169 0.258819075248021502 0 1 6.12323399573676604e-17 0 1.22464679914735321e-16 -1 0 -0.25881916568452229 -0.965925793979210301 0 -0.499999891889129133 -0.866025466202270056 0 -0.707106781186547462 -0.707106781186547684 0 -0.866025466202269723 -0.499999891889129799 0 -0.965925793979210301 -0.258819165684522123 0 -1 -1.83697019872102969e-16 0 -0.965925793979210412 0.25881916568452179 0 -0.866025466202269945 0.499999891889129411 0 -0.707106781186547684 0.707106781186547462 0 -0.499999891889129799 0.866025466202269723 0 -0.258819165684522179 0.965925793979210301 0 0.5 3.06161699786838302e-17 0 0.25 1.53080849893419151e-17 0 0.75 4.59242549680257422e-17 0 6.12323399573676604e-17 -0.5 0 3.06161699786838302e-17 -0.25 0 9.18485099360514844e-17 -0.75 0 -0.366742427287351624 0.49448770056279745 0 -0.543519020850925627 -0.0930502013393769284 0 0.271588156561690552 0.571569303201864054 0 -0.333187186624814335 -0.51024527210371351 0 -0.673915622739848508 0.283561802523850326 0 0.607004670678171987 0.391362141271555775 0 -0.0794826655388244929 0.728066213011798524 0 -0.23631871621712236 -0.141029353317234207 0 -0.297582245224840836 0.137725505862699771 0 -0.733530383789228524 -0.340506243303687339 0 0.0228642765842100235 0.300291579746488013 0 0.332678628136344745 0.244418372194045025 0</vertices>
+	<edges>2 1 3 2 4 3 5 4 6 5 7 6 9 8 10 9 11 10 12 11 13 12 14 13 15 14 16 15 17 16 18 17 19 18 1 19 0 21 21 20 20 22 22 7 0 24 24 23 23 25 25 8 22 6 9 25 17 26 26 18 27 14 2 28 3 28 23 29 29 11 27 29 10 29 24 29 25 29 9 29 15 30 30 26 27 30 16 30 14 30 17 30 5 31 31 20 28 31 22 31 4 31 6 31 3 31 26 32 32 1 28 32 19 32 2 32 18 32 0 33 33 24 27 33 29 33 27 34 34 26 33 34 30 34 0 34 12 35 35 27 29 35 13 35 11 35 14 35 21 36 36 0 28 36 26 36 34 36 32 36 31 37 37 20 36 37 28 37 21 37</edges>
+	<triangles>7 6 22 8 25 9 18 17 26 3 2 28 11 10 29 23 24 29 25 23 29 9 25 29 10 9 29 16 15 30 14 27 30 15 14 30 17 16 30 26 17 30 20 22 31 5 4 31 6 5 31 22 6 31 3 28 31 4 3 31 1 19 32 2 1 32 28 2 32 18 26 32 19 18 32 24 0 33 29 24 33 27 29 33 27 33 34 30 27 34 26 30 34 33 0 34 29 27 35 13 12 35 11 29 35 12 11 35 14 13 35 27 14 35 0 21 36 34 0 36 26 34 36 32 26 36 28 32 36 20 31 37 28 36 37 31 28 37 21 20 37 36 21 37</triangles>
+	<subset_handler name="defSH">
+		<subset name="Inner" color="0.588235 0.588235 1 1" state="262144">
+			<vertices>26 27 28 29 30 31 32 33 34 35 36 37</vertices>
+			<edges>26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84</edges>
+			<faces>0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47</faces>
+		</subset>
+		<subset name="Circle" color="1 0 0 1" state="262144">
+			<vertices>1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19</vertices>
+			<edges>0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17</edges>
+		</subset>
+		<subset name="Cut" color="0 1 0 1" state="0">
+			<vertices>0 20 21 22 23 24 25</vertices>
+			<edges>18 19 20 21 22 23 24 25</edges>
+		</subset>
+	</subset_handler>
+	<subset_handler name="markSH">
+		<subset name="crease" color="1 1 1 1" state="0"/>
+		<subset name="fixed" color="1 1 1 1" state="0"/>
+	</subset_handler>
+	<selector name="defSel">
+		<vertices>21 1</vertices>
+	</selector>
+	<projection_handler name="defPH" subset_handler="0">
+		<default type="default">0 0</default>
+		<projector type="sphere" subset="1">0 0 0.00000000000000000e+00 0.00000000000000000e+00 0.00000000000000000e+00 -1.00000000000000000e+00 -1.00000000000000000e+00</projector>
+	</projection_handler>
+</grid>
+
diff --git a/content/tutorial-fem/grids/unit_square_quad.ugx b/content/tutorial-fem/grids/unit_square_quad.ugx
new file mode 100644
index 0000000..0f8db9a
--- /dev/null
+++ b/content/tutorial-fem/grids/unit_square_quad.ugx
@@ -0,0 +1,18 @@
+<?xml version="1.0" encoding="utf-8"?>
+<grid name="defGrid">
+	<vertices coords="3">1 1 0 1 0 0 0 0 0 0 1 0 0.5 1 0 0 0.5 0 0.5 0 0 1 0.5 0 0.5 0.5 0</vertices>
+	<edges>0 4 4 3 3 5 5 2 2 6 6 1 1 7 7 0 4 8 8 7 5 8 6 8</edges>
+	<quadrilaterals>0 4 8 7 3 5 8 4 2 6 8 5 1 7 8 6</quadrilaterals>
+	<subset_handler name="defSH">
+		<subset name="Inner" color="1 0 0 1">
+			<vertices>8</vertices>
+			<edges>8 9 10 11</edges>
+			<faces>0 1 2 3</faces>
+		</subset>
+		<subset name="Boundary" color="0 1 0 1">
+			<vertices>0 3 2 1 4 5 6 7</vertices>
+			<edges>0 1 2 3 4 5 6 7</edges>
+		</subset>
+	</subset_handler>
+</grid>
+
diff --git a/content/tutorial-fem/grids/unit_square_tri.ugx b/content/tutorial-fem/grids/unit_square_tri.ugx
new file mode 100644
index 0000000..4fac420
--- /dev/null
+++ b/content/tutorial-fem/grids/unit_square_tri.ugx
@@ -0,0 +1,29 @@
+<?xml version="1.0" encoding="utf-8"?>
+<grid name="defGrid">
+	<vertices coords="3">1 1 0 1 0 0 0 0 0 0 1 0 0.5 1 0 0 0.5 0 0.5 0 0 1 0.5 0 0.5 0.5 0</vertices>
+	<edges>0 4 4 3 3 5 5 2 2 6 6 1 1 7 7 0 4 8 8 7 5 8 6 8 8 0 8 2 5 4 7 6</edges>
+	<triangles>0 4 8 8 7 0 2 6 8 8 5 2 5 8 4 3 5 4 7 8 6 1 7 6</triangles>
+	<subset_handler name="defSH">
+		<subset name="Inner" color="1 0 0 1" state="262144">
+			<vertices>8</vertices>
+			<edges>8 9 10 11 12 13 14 15</edges>
+			<faces>0 1 2 3 4 5 6 7</faces>
+		</subset>
+		<subset name="Boundary" color="0 1 0 1" state="393216">
+			<vertices>0 3 2 1 4 5 6 7</vertices>
+			<edges>0 1 2 3 4 5 6 7</edges>
+		</subset>
+	</subset_handler>
+	<subset_handler name="markSH">
+		<subset name="crease" color="1 1 1 1" state="0"/>
+		<subset name="fixed" color="1 1 1 1" state="0"/>
+	</subset_handler>
+	<selector name="defSel">
+		<edges>15 1 1 1</edges>
+		<faces>5 1 4 1 0 1 1 1 6 1 7 1 2 1 3 1</faces>
+	</selector>
+	<projection_handler name="defPH" subset_handler="0">
+		<default type="default">0 0</default>
+	</projection_handler>
+</grid>
+