diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
index 9ed876922d76c25872f7f918500c77bbb6893c4f..5cfa0fb5d79bd257bcab583a98d0052435d4e29d 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.cpp
@@ -23,6 +23,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
         const BoundaryConditionConfig& config,
         const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
         const int variable_id, const unsigned integration_order,
+        const unsigned shapefunction_order,
         const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
         MeshGeoToolsLib::MeshNodeSearcher& /*mesh_node_searcher*/,
         MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher)
@@ -31,7 +32,7 @@ BoundaryConditionBuilder::createNeumannBoundaryCondition(
         config.config,
         getClonedElements(boundary_element_searcher, config.geometry),
         dof_table, variable_id, config.component_id,
-        mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+        mesh.isAxiallySymmetric(), integration_order, shapefunction_order, mesh.getDimension(),
         parameters, _fracture_prop);
 }
 
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
index 20fe4f15fde80725247604c25c87f013c653c91f..5e7bac36b86e2bce00a1ef14850703f6f6c24286 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/BoundaryConditionBuilder.h
@@ -45,6 +45,7 @@ private:
         const NumLib::LocalToGlobalIndexMap& dof_table,
         const MeshLib::Mesh& mesh, const int variable_id,
         const unsigned integration_order,
+        const unsigned shapefunction_order,
         const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
             parameters,
         MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher,
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 2eaa6b10edb31225a8924c33ff41b9dad683244b..84222f6cd0b8a0b125fd249fe5eaa0b8627f122d 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -33,6 +33,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
                          typename std::decay<Data>::type>::value,
             bool>::type is_axially_symmetric,
         unsigned const integration_order,
+        unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
         unsigned const global_dim,
@@ -63,7 +64,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
         variable_id, component_id, std::move(all_mesh_subsets), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
-        global_dim, _elements, *_dof_table_boundary, _local_assemblers,
+        global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
         is_axially_symmetric, _integration_order, _data, fracture_prop, variable_id);
 }
 
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 3144eb2b31520aa4cc24cf8e69b8f8b53b772a1b..4a2e677eb794b698a5cbd9fd82c7a0ba3fb5e63b 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -37,6 +37,7 @@ public:
                          typename std::decay<Data>::type>::value,
             bool>::type is_axially_symmetric,
         unsigned const integration_order,
+        unsigned const shapefunction_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
         unsigned const global_dim,
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
index 5f99aedf0147c6b104958a94b005f8eb3767e50f..fb09e7c3538fd4e9bd1c6aa99482ccebc3ef57f6 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -22,13 +22,13 @@ namespace SmallDeformationWithLIE
 using NeumannBoundaryCondition = GenericNaturalBoundaryCondition<
     Parameter<double> const&, NeumannBoundaryConditionLocalAssembler>;
 
-std::unique_ptr<BoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<BoundaryCondition> createNeumannBoundaryCondition(
     BaseLib::ConfigTree const& config,
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-        int const component_id, bool is_axially_symmetric,
-        unsigned const integration_order, unsigned const global_dim,
+    int const component_id, bool is_axially_symmetric,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     FractureProperty const& fracture_prop)
 {
@@ -44,7 +44,7 @@ createNeumannBoundaryCondition(
 
     return std::unique_ptr<BoundaryCondition>(
         new NeumannBoundaryCondition(
-            is_axially_symmetric, integration_order,
+            is_axially_symmetric, integration_order, shapefunction_order,
             dof_table, variable_id, component_id,
             global_dim, std::move(elements), param, fracture_prop));
 }
diff --git a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
index 5c846581e8b89df312c288e9e2bbf1a8a1f74d41..dc76409af349e8d4a905d7c52f0da3a282eff17f 100644
--- a/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/SmallDeformationWithLIE/BoundaryCondition/NeumannBoundaryCondition.h
@@ -32,7 +32,9 @@ createNeumannBoundaryCondition(
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
     int const component_id, bool is_axially_symmetric,
-    unsigned const integration_order, unsigned const global_dim,
+    unsigned const integration_order,
+    unsigned const shapefunction_order,
+    unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     FractureProperty const& fracture_prop);