diff --git a/tests/3d_r13/inputs/r13_kn0.005_shell.yml b/tests/3d_r13/inputs/r13_kn0.005_shell.yml
index 2fc35b14bdcfd8c10615a31305d9c0b7bbe13254..b6bd902de525ab14c6e29d35ef2423da2ade01b9 100644
--- a/tests/3d_r13/inputs/r13_kn0.005_shell.yml
+++ b/tests/3d_r13/inputs/r13_kn0.005_shell.yml
@@ -10,8 +10,8 @@ meshes:
   - ../3d_mesh/shell1.h5
   - ../3d_mesh/shell2.h5
   # - ../3d_mesh/shell3.h5
-  #- ../3d_mesh/shell4.h5
-  #- ../3d_mesh/shell5.h5
+  # - ../3d_mesh/shell4.h5
+  # - ../3d_mesh/shell5.h5
 
 # Numerical Parameters
 # ====================
diff --git a/tests/3d_r13/inputs/r13_kn0.2_shell.yml b/tests/3d_r13/inputs/r13_kn0.2_shell.yml
new file mode 100644
index 0000000000000000000000000000000000000000..518b0cd7b284bcbac5491197b989332515696505
--- /dev/null
+++ b/tests/3d_r13/inputs/r13_kn0.2_shell.yml
@@ -0,0 +1,196 @@
+# General
+# =======
+# - output_folder: Used as output folder
+output_folder: r13_kn0.2_shell
+
+# Meshes
+# ======
+# - meshes: List of input meshes in h5 format to run simulations on
+meshes:
+  - ../3d_mesh/shell1.h5
+  - ../3d_mesh/shell2.h5
+  # - ../3d_mesh/shell3.h5
+  # - ../3d_mesh/shell4.h5
+  # - ../3d_mesh/shell5.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 heatflu_x_w 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: False
+    delta_theta: 1
+    delta_u: 1
+    delta_p: 1
+  gls:
+    enable: True
+    tau_energy: 0.1
+    tau_heatflux: 0.1
+    tau_mass: 0.01
+    tau_momentum: 10
+    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
+# - f_s: Heatflux force for mode==heat||r13
+# - f_sigma: Stress force for mode==stress||r13
+nsd: 3
+mode: r13
+heat_source: 0
+mass_source: 0
+body_force: [0,0,0]
+f_s: [0,0,0]
+f_sigma: [[0,0,0],[0,0,0],[0,0,0]]
+
+# PETSc Options (dictionary)
+# ==========================
+# [#1 e.g. for mumps directy solver:]
+# - ksp_type: preonly  # preconditioner only, i.e. 1 application
+# - pc_type: lu  # use LU as preconditioner <=> direct solve
+# - pc_factor_mat_solver_type: mumps  # use mumps sparse solver
+# [#2 e.g. for gmres iterative solver with icc preconditioner:]
+# - ksp_type: gmres  # Generalized Minimal Residual Method
+# - pc_type: icc  # incomplete Cholesky
+petsc_options:
+  ksp_type: tfqmr
+  ksp_rtol: 1E-7
+  ksp_max_it: 100000
+  ksp_view:
+  ksp_monitor_true_residual:
+  pc_type: fieldsplit
+  pc_fieldsplit_detect_saddle_point:
+  # option 1:
+  pc_fieldsplit_type: schur
+  pc_fieldsplit_schur_fact_type: full
+  pc_fieldsplit_schur_precondition: selfp
+  fieldsplit_0_ksp_type: preonly
+  fieldsplit_0_pc_type: gamg
+  fieldsplit_0_pc_gamg_type: classical
+  fieldsplit_1_ksp_type: preonly
+  fieldsplit_1_pc_type: jacobi
+
+# Region Parameters
+# =================
+# - regs: Dictionary of all mesh regions
+#   - reg_id: Must contain the following parameters:
+#     - kn: Knudsen number
+regs:
+  10002: #inner mesh region
+    kn: 0.2
+
+# Boundary Conditions
+# ===================
+# - polar_coord_syst: true needs u_n_w,u_t_w; false needs u_x_w,u_y_w,u_z_w
+# - 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 (for polar_coord_syst=true)
+#     - u_n_w: Value for normal velocity at wall (for polar_coord_syst=true)
+#     - u_x_w: Value for x-velocity at wall (for polar_coord_syst=false)
+#     - u_y_w: Value for y-velocity at wall (for polar_coord_syst=false)
+#     - u_z_w: Value for z-velocity at wall (for polar_coord_syst=false&&nsd=3)
+#     - p_w: Value for pressure at wall
+#     - epsilon_w: Inflow-model parameter <=> Weight of pressure prescription
+polar_coord_syst: False
+bcs:
+  10000: # inner
+    chi_tilde: 1.0
+    theta_w: 1
+    u_t_w: 1E300
+    u_n_w: 1E300
+    u_x_w: 0
+    u_y_w: 0
+    u_z_w: 0
+    p_w: 0
+    epsilon_w: 0.01
+  10001: # outer
+    chi_tilde: 1.0
+    theta_w: 2
+    u_t_w: 1E300
+    u_n_w: 1E300
+    u_x_w: 0
+    u_y_w: 0
+    u_z_w: 1
+    p_w: -0.027*x[2]/pow((pow(x[0],2)+pow(x[1],2)+pow(x[2],2)),0.5)
+    epsilon_w: 100
+
+# 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: Shift numerical pressure (False,zeromean,zerominimum)
+# - relative_errors: Use relative errors. If exact sol. is zero, use absolute.
+convergence_study:
+  enable: True
+  exact_solution: esols/R13_Kn0.2_eps0.01_p0.027.cpp
+  plot: False # to avoid error exit code due to $DISPLAY
+  write_systemmatrix: False
+  rescale_pressure: zeromean
+  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: []
+
+# 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: []