Added functionality to specify degree element order

src/festim/hydrogen_transport_problem.py

@@ -274,7 +274,7 @@ def volume_meshtags(self, value):
 
     def initialise(self):
         self.create_species_from_traps()
-        self.define_function_spaces()
+        self.define_function_spaces(element_degree=self.settings.element_degree)
         self.define_meshtags_and_measures()
         self.assign_functions_to_species()
 
@@ -531,25 +531,27 @@ def define_D_global(self, species):
         D.interpolate(D_expr)
         return D, D_expr
 
-    def define_function_spaces(self):
+    def define_function_spaces(self, element_degree=1):
         """Creates the function space of the model, creates a mixed element if
         model is multispecies. Creates the main solution and previous solution
         function u and u_n. Create global DG function spaces of degree 0 and 1
-        for the global diffusion coefficient"""
+        for the global diffusion coefficient.
+
+        Args:
+            element_degree (int, optional): Degree order for finite element.
+                Defaults to 1.
+        """
 
-        # TODO: expose degree as a property to the user (element_degree ?)
-        # in ProblemBase
-        degree = 1
         element_CG = basix.ufl.element(
             basix.ElementFamily.P,
             self.mesh.mesh.basix_cell(),
-            degree,
+            element_degree,
             basix.LagrangeVariant.equispaced,
         )
         element_DG = basix.ufl.element(
             "DG",
             self.mesh.mesh.basix_cell(),
-            degree,
+            element_degree,
             basix.LagrangeVariant.equispaced,
         )
 
@@ -1143,7 +1145,9 @@ def create_initial_conditions(self):
                 "initial conditions not yet implemented for discontinuous"
             )
 
-    def define_function_spaces(self, subdomain: _subdomain.VolumeSubdomain):
+    def define_function_spaces(
+        self, subdomain: _subdomain.VolumeSubdomain, element_degree=1
+    ):
         """
         Creates appropriate function space and functions for a given subdomain (submesh)
         based on the number of species existing in this subdomain. Then stores the
@@ -1154,6 +1158,8 @@ def define_function_spaces(self, subdomain: _subdomain.VolumeSubdomain):
 
         Args:
             subdomain (F.VolumeSubdomain): a subdomain of the geometry
+            element_degree (int, optional): Degree order for finite element.
+                Defaults to 1.
         """
         # get number of species defined in the subdomain
         all_species = [
@@ -1167,11 +1173,10 @@ def define_function_spaces(self, subdomain: _subdomain.VolumeSubdomain):
                 unique_species.append(species)
         nb_species = len(unique_species)
 
-        degree = 1
         element_CG = basix.ufl.element(
             basix.ElementFamily.P,
             subdomain.submesh.basix_cell(),
-            degree,
+            element_degree,
             basix.LagrangeVariant.equispaced,
         )
         element = basix.ufl.mixed_element([element_CG] * nb_species)

src/festim/settings.py

@@ -14,6 +14,8 @@ class Settings:
         transient (bool, optional): Whether the simulation is transient or not.
         final_time (float, optional): Final time for a transient simulation.
             Defaults to None
+        element_degree (int, optional): Degree order for finite element.
+            Defaults to 1.
         stepsize (festim.Stepsize, optional): stepsize for a transient
             simulation. Defaults to None
         convergence_criterion: resiudal or incremental (for Newton solver)
@@ -24,6 +26,7 @@ class Settings:
         max_iterations (int): Maximum number of iterations for the solver.
         transient (bool): Whether the simulation is transient or not.
         final_time (float): Final time for a transient simulation.
+        element_degree (int): Degree order for finite element.
         stepsize (festim.Stepsize): stepsize for a transient
             simulation.
         convergence_criterion: resiudal or incremental (for Newton solver)
@@ -37,6 +40,7 @@ def __init__(
         max_iterations=30,
         transient=True,
         final_time=None,
+        element_degree=1,
         stepsize=None,
         convergence_criterion: Literal["residual", "incremental"] = "residual",
     ) -> None:
@@ -45,6 +49,7 @@ def __init__(
         self.max_iterations = max_iterations
         self.transient = transient
         self.final_time = final_time
+        self.element_degree = element_degree
         self.stepsize = stepsize
         self.convergence_criterion = convergence_criterion