diff --git a/examples/thermal_edge_flow/create_mesh9.sh b/examples/thermal_edge_flow/create_mesh9.sh
new file mode 100755
index 0000000000000000000000000000000000000000..6a0934e653d36c5738bf060ad1942ceabba14a15
--- /dev/null
+++ b/examples/thermal_edge_flow/create_mesh9.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+
+for split in $(seq 0 1 3)
+do
+  outname=study9_"$split".h5
+  echo $outname
+  geoToH5 study9.geo "$outname" "-setnumber split $split"
+done
diff --git a/examples/thermal_edge_flow/input9.yml b/examples/thermal_edge_flow/input9.yml
new file mode 100644
index 0000000000000000000000000000000000000000..f903175662c5ec0e12a6d560b258e5438941a668
--- /dev/null
+++ b/examples/thermal_edge_flow/input9.yml
@@ -0,0 +1,165 @@
+# General
+# =======
+# - output_folder: Used as output folder
+# output_folder:
+output_folder: study9
+
+# Meshes
+# ======
+# - meshes: List of input meshes in h5 format to run simulations on
+meshes:
+  - study9_0.h5
+  - study9_1.h5
+  - study9_2.h5
+  - study9_3.h5
+
+# Numerical Parameters
+# ====================
+# - elements: Must contain the fields: theta, s, p, u, sigma
+#   - fields: List of FEM parameters (shape, degree)
+#     - shape: Element shape, e.g. Lagrange
+#     - degree: Element degree, e.g. 2
+# - stabilization: Must contain cip and gls
+#   - cip: Collection of Continous Interior Penalty (CIP) parameters
+#     - enable: Enable CIP stabilization
+#     - delta_theta: Stabilization of grad(T)*grad(T_test) over edge
+#     - delta_u: Stabilization of grad(u)*grad(u_test) over edge
+#     - delta_p: Stabilization of grad(p)*grad(p_test) over edge
+#   - gls: Collection of Garlerkin Least Squares (GLS) parameters
+#     - enable: Enable GLS stabilization
+#     - tau_energy: Stabilization with energy eq. residual
+#     - tau_heatflux: Stabilization with heatflux eq. residual
+#     - tau_mass: Stabilization with mass eq. residual
+#     - tau_momentum: Stabilization with momentum eq. residual
+#     - tau_stress: Stabilization with stress eq. residual
+elements:
+  theta:
+    shape: Lagrange
+    degree: 1
+  s:
+    shape: Lagrange
+    degree: 1
+  p:
+    shape: Lagrange
+    degree: 1
+  u:
+    shape: Lagrange
+    degree: 1
+  sigma:
+    shape: Lagrange
+    degree: 1
+stabilization:
+  cip:
+    enable: True
+    delta_theta: 1.0
+    delta_u: 1.0
+    delta_p: 0.1
+  gls:
+    enable: False
+    tau_energy: 0.001
+    tau_heatflux: 0.001
+    tau_mass: 0.01
+    tau_momentum: 0.01
+    tau_stress: 0.01
+
+# Formulation Parameters
+# ======================
+# - nsd: Number of spatial dimensions == 2
+# - mode: Formulation mode, one of heat, stress, r13
+# - heat_source: Heat source function for mode==heat||r13
+# - mass_source: Mass source function for mode==stress||r13
+# - body_force: Body force for mode==stress||r13
+nsd: 2
+mode: r13
+heat_source: 0
+mass_source: 0
+body_force: [0,0]
+
+# Region Parameters
+# =================
+# - regs: Dictionary of all mesh regions
+#   - reg_id: Must contain the following parameters:
+#     - kn: Knudsen number
+regs:
+  4000:
+    # kn: 0.003544908 # 2*1.772453851*0.001
+    kn: 0.001
+
+# Boundary Conditions
+# ===================
+# - bcs: Dictionary of all boundary IDs from mesh
+#   - bc_id: must contain the following parameters
+#     - chi_tilde: Refaction coefficient in Maxwell accomodation model
+#     - theta_w: Value for temperature at wall
+#     - u_t_w: Value for tangential velocity at wall
+#     - u_n_w: Value for normal velocity at wall
+#     - p_w: Value for pressure at wall
+#     - epsilon_w: Inflow-model parameter <=> Weight of pressure prescription
+bcs:
+  3000: # outer
+    chi_tilde: 1.0
+    theta_w: 0
+    u_t_w: 0
+    u_n_w: 0
+    p_w: 0
+    epsilon_w: 0
+  3100: # inner
+    chi_tilde: 1.0
+    theta_w: 1
+    u_t_w: 0
+    u_n_w: 0
+    p_w: 0
+    epsilon_w: 0
+
+# Convergence Study
+# =================
+# - enable: Enable convergence study on given meshes
+# - exact_solution: Path to exact solution in cpp-format to compare errors
+# - plot: Show errors in matplotlib window. PDF output is always per default.
+# - write_systemmatrix: Writes out systemmatrix (LHS) to use for analysis
+# - rescale_pressure: Rescale numerical pressure result to have zero mean
+# - relative_errors: Use relative errors. If exact sol. is zero, use absolute.
+convergence_study:
+  enable: False
+  exact_solution: esols/01_coeffs.cpp
+  plot: False # to avoid error exit code due to $DISPLAY
+  write_systemmatrix: False
+  rescale_pressure: True
+  relative_error: True
+
+# Postprocessing
+# ==============
+# - write_pdfs: Write all solution fields as PDF plot
+# - write_vecs: Write all solution fields as vectors
+# - massflow: List of BC IDs to compute massflow J=int_bc dot(u,n) ds
+# - line_integrals: List of line integral dicts:
+#   - name: Name for output
+#   - expr: Expression to evaluate
+#   - start: Start point
+#   - end: End point
+#   - res: Sampling resolution of line
+postprocessing:
+  write_pdfs: False
+  write_vecs: False
+  massflow: []
+  line_integrals:
+    - name: "avg(abs(uy(x,y=0.5)))"
+      expr: abs(uy)/8
+      start: [0, 0.5]
+      end: [8, 0.5]
+      res: 10000
+    - name: "avg(abs(uy(x,y=4.5)))"
+      expr: abs(uy)/8
+      start: [0, 4.5]
+      end: [8, 4.5]
+      res: 10000
+
+# Parameter Study
+# ==============
+# - enable: Repeat simulation with different p. values (study)
+# - parameter_key: Key as list, e.g. ["elemenets", "p", "degree"]
+# - parameter_values: List of value for parameter, e.g. [0.01,0.1,1,10]
+parameter_study:
+  enable: False
+  parameter_key: []
+  parameter_values: []
diff --git a/examples/thermal_edge_flow/slurm.sh b/examples/thermal_edge_flow/slurm.sh
new file mode 100644
index 0000000000000000000000000000000000000000..33d882ffbebd7acdea0dda6ab8da635f07f26c69
--- /dev/null
+++ b/examples/thermal_edge_flow/slurm.sh
@@ -0,0 +1,10 @@
+#!/usr/local_rwth/bin/zsh
+
+#SBATCH --ntasks=100
+#SBATCH --time=0-01:00:00
+#SBATCH --mem-per-cpu=30000M
+#SBATCH --job-name=fenicsR13_test
+#SBATCH --output=output.%J.txt
+
+cd $HOME/fenicsR13/examples/thermal_edge_flow
+$MPIEXEC $FLAGS_MPI_BATCH singularity exec /rwthfs/rz/SW/UTIL.common/singularity/fenicsr13_latest.sif fenicsR13 input9.yml
diff --git a/examples/thermal_edge_flow/study9.geo b/examples/thermal_edge_flow/study9.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3a5aa0d60fc7fcec474ec77ab4a504470d497acb
--- /dev/null
+++ b/examples/thermal_edge_flow/study9.geo
@@ -0,0 +1,147 @@
+Mesh.MshFileVersion = 2.0;
+
+// Command line Parameters
+
+If(!Exists(split))
+  split = 0;
+EndIf
+Printf("split=%g", split);
+
+
+dist1 = 0.12;
+dist2 = 0.4;
+
+res_tmp1 = 2^-split * 2^-6; // leftbot edge
+res_tmp2 = 2^-split * 2^-3; // bulk
+res_tmp3 = 2^-split * 2^-4; // topright edge
+res_tmp4 = 2^-split * 2^-7; // leftbot focus
+res_tmp5 = 2^-split * 2^-11; // inner edge focus
+res_tmp6 = 2^-split * 2^-7; // inner edge
+res_tmp7 = 2^-split * 2^-4; // focus bulk
+
+
+Point(1001) = {0.0, 0.0, 0, res_tmp1};
+Point(1002) = {2.0, 0.0, 0, res_tmp4};
+Point(1003) = {4.0, 0.0, 0, res_tmp1};
+Point(1004) = {8.0-dist1, 0.0, 0, res_tmp3};
+Point(1005) = {8.0, 0.0, 0, res_tmp3};
+Point(1006) = {8.0, 0.0+dist1, 0, res_tmp3};
+Point(1007) = {8.0, 8.0-dist1, 0, res_tmp3};
+Point(1008) = {8.0, 8.0, 0, res_tmp3};
+Point(1009) = {8.0-dist1, 8.0, 0, res_tmp3};
+Point(1010) = {0.0+dist1, 8.0, 0, res_tmp3};
+Point(1011) = {0.0, 8.0, 0, res_tmp3};
+Point(1012) = {0.0, 8.0-dist1, 0, res_tmp3};
+Point(1013) = {0.0, 4.0, 0, res_tmp1};
+Point(1014) = {0.0, 2.0, 0, res_tmp4};
+
+Point(1101) = {0.0+dist1, 0.0+dist1, 0, res_tmp7};
+Point(1102) = {4.0, 0.0+dist1, 0, res_tmp7};
+Point(1103) = {8.0-dist1, 0.0+dist1, 0, res_tmp2};
+Point(1104) = {8.0-dist1, 8.0-dist1, 0, res_tmp2};
+Point(1105) = {0.0+dist1, 8.0-dist1, 0, res_tmp2};
+Point(1106) = {0.0+dist1, 4.0, 0, res_tmp7};
+
+Point(1202) = {1.0, 1.0-dist2, 0, res_tmp7};
+Point(1203) = {3.0, 1.0-dist2, 0, res_tmp7};
+Point(1205) = {3.0+dist2, 1.0, 0, res_tmp7};
+Point(1206) = {3.0+dist2, 3.0, 0, res_tmp7};
+Point(1208) = {3.0, 3.0+dist2, 0, res_tmp7};
+Point(1209) = {1.0, 3.0+dist2, 0, res_tmp7};
+Point(1211) = {1.0-dist2, 3.0, 0, res_tmp7};
+Point(1212) = {1.0-dist2, 1.0, 0, res_tmp7};
+
+Point(1301) = {1.0, 1.0, 0, res_tmp5};
+Point(1302) = {1.0+dist2, 1.0, 0, res_tmp6};
+Point(1303) = {3.0-dist2, 1.0, 0, res_tmp6};
+Point(1304) = {3.0, 1.0, 0, res_tmp5};
+Point(1305) = {3.0, 1.0+dist2, 0, res_tmp6};
+Point(1306) = {3.0, 3.0-dist2, 0, res_tmp6};
+Point(1307) = {3.0, 3.0, 0, res_tmp5};
+Point(1308) = {3.0-dist2, 3.0, 0, res_tmp6};
+Point(1309) = {1.0+dist2, 3.0, 0, res_tmp6};
+Point(1310) = {1.0, 3.0, 0, res_tmp5};
+Point(1311) = {1.0, 3.0-dist2, 0, res_tmp6};
+Point(1312) = {1.0, 1.0+dist2, 0, res_tmp6};
+
+
+Point(1401) = {4.0, 4.0, 0, res_tmp7};
+
+
+Line(2001) = {1001, 1002};
+Line(2002) = {1002, 1003};
+Line(2003) = {1003, 1004};
+Line(2004) = {1004, 1005};
+Line(2005) = {1005, 1006};
+Line(2006) = {1006, 1007};
+Line(2007) = {1007, 1008};
+Line(2008) = {1008, 1009};
+Line(2009) = {1009, 1010};
+Line(2010) = {1010, 1011};
+Line(2011) = {1011, 1012};
+Line(2012) = {1012, 1013};
+Line(2013) = {1013, 1014};
+Line(2014) = {1014, 1001};
+
+
+Line(2101) = {1101, 1102};
+Line(2102) = {1102, 1103};
+Line(2103) = {1103, 1104};
+Line(2104) = {1104, 1105};
+Line(2105) = {1105, 1106};
+Line(2106) = {1106, 1101};
+
+
+Line(2202) = {1202, 1203};
+Line(2205) = {1205, 1206};
+Line(2208) = {1208, 1209};
+Line(2211) = {1211, 1212};
+
+Line(2301) = {1301, 1302};
+Line(2302) = {1302, 1303};
+Line(2303) = {1303, 1304};
+Line(2304) = {1304, 1305};
+Line(2305) = {1305, 1306};
+Line(2306) = {1306, 1307};
+Line(2307) = {1307, 1308};
+Line(2308) = {1308, 1309};
+Line(2309) = {1309, 1310};
+Line(2310) = {1310, 1311};
+Line(2311) = {1311, 1312};
+Line(2312) = {1312, 1301};
+
+Line(2401) = {1106, 1401};
+Line(2402) = {1401, 1102};
+
+Circle(2601) = {1212, 1301, 1202};
+Circle(2602) = {1203, 1304, 1205};
+Circle(2603) = {1206, 1307, 1208};
+Circle(2604) = {1209, 1310, 1211};
+
+
+Curve Loop(3207) = {2202, 2602, 2205, 2603, 2208, 2604, 2211, 2601};
+Curve Loop(3208) = {2301:2312};
+Plane Surface(4207) = {3207, 3208};
+Curve Loop(3209) = {2105, 2106, 2101, 2102, 2103, 2104};
+Curve Loop(3210) = {2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2001};
+Plane Surface(4209) = {3210, 3209};
+Curve Loop(3211) = {2102, 2103, 2104, 2105, 2401, 2402};
+Plane Surface(4210) = {3211};
+Curve Loop(3212) = {2101, -2402, -2401, 2106};
+Plane Surface(4211) = {3212, 3207};
+
+
+Physical Curve("outer",3000) = {
+  2001:2014
+};
+
+Physical Curve("inner",3100) = {
+  2301:2312
+};
+
+Physical Surface("mesh",4000) = {
+  4207,
+  4209,
+  4210,
+  4211
+};