diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index baba8990961f9f40480a4b4ae35c710413c1300b..9959bbfcd0d21e15ef3f9aa751f55a5a887458f6 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -33,7 +33,8 @@ class LDDsimulation(object):
                  tol: float = None,#
                  LDDsolver_tol: float = 1E-8,
                  # maximal number of L-iterations that the LDD solver uses.
-                 max_iter_num: int = 1000):
+                 max_iter_num: int = 1000,
+                 FEM_Lagrange_degree: int = 1):
         hlp.print_once("\n ###############################################\n",
                        "############# Begin Simulation ################\n",
                        "###############################################\n")
@@ -48,6 +49,7 @@ class LDDsimulation(object):
         # dimension of the simulation. This is by default 2 as I don't intend
         # to implement 3D.
         self.dim = dimension
+        self.FEM_Lagrange_degree = FEM_Lagrange_degree
         #Parameters
         # df.parameters["allow_extrapolation"] = True
         # df.parameters["refinement_algorithm"] = "plaza_with_parent_facets"
@@ -1107,7 +1109,8 @@ class LDDsimulation(object):
                         saturation = self.saturation[subdom_num],#
                         timestep_size = self.timestep_size,#
                         interpolation_degree = self.interpolation_degree,#
-                        tol = self.tol)
+                        tol = self.tol,
+                        degree=self.FEM_Lagrange_degree)
                     }
                 )
             # setup the linear solvers and set the solver parameters.