diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 049584c707909560d0f9b9e539b8c7c55b919e8a..ece058652cab9c7b95924506e5534f02731bc08b 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -22,6 +22,7 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
     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 =
@@ -37,21 +38,21 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
     {
         return createDirichletBoundaryCondition(
                     config, dof_table, mesh, variable_id,
-                    integration_order, parameters,
+                    integration_order, shapefunction_order, parameters,
                     mesh_node_searcher, boundary_element_searcher);
     }
     else if (type == "Neumann")
     {
         return createNeumannBoundaryCondition(
                     config, dof_table, mesh, variable_id,
-                    integration_order, parameters,
+                    integration_order, shapefunction_order, parameters,
                     mesh_node_searcher, boundary_element_searcher);
     }
     else if (type == "Robin")
     {
         return createRobinBoundaryCondition(
                     config, dof_table, mesh, variable_id,
-                    integration_order, parameters,
+                    integration_order, shapefunction_order, parameters,
                     mesh_node_searcher, boundary_element_searcher);
     }
     else
@@ -65,6 +66,7 @@ BoundaryConditionBuilder::createDirichletBoundaryCondition(
         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*/)
@@ -110,6 +112,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)
@@ -118,7 +121,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);
 }
 
@@ -127,6 +130,7 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
         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)
@@ -135,7 +139,7 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
         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);
 }
 
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 32f8df74df5e801eb6e9a277a6cd15e95651cdb8..3a6b78039059f82fb8f0a0404423344501bf0cda 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -76,6 +76,7 @@ public:
         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);
 
 protected:
@@ -84,6 +85,7 @@ protected:
         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,
@@ -94,6 +96,7 @@ protected:
         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,
@@ -104,6 +107,7 @@ protected:
         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/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 515b653b8592c1332358f0e91c71e93cc1e2a89b..25e8988090a56fc94d5194228ff59f66ab5508a3 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -26,6 +26,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, std::vector<MeshLib::Element*>&& elements,
@@ -55,7 +56,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);
 }
 
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index 5edcf191a7302fc8a42439aa265a22fa8700a902..601ec399fcd15c70e8c58ce66d2ac7a141cac3c5 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -33,6 +33,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, std::vector<MeshLib::Element*>&& elements,
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
index 79601405d89f4dc9446a89cdb636d37a45a83a74..7a4f13dfe295afcb166327c7c9e1d12afe519d00 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -17,7 +17,9 @@ std::unique_ptr<NeumannBoundaryCondition> 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)
 {
     DBUG("Constructing Neumann BC from config.");
@@ -31,7 +33,7 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     auto const& param = findParameter<double>(param_name, parameters, 1);
 
     return std::unique_ptr<NeumannBoundaryCondition>(
-        new NeumannBoundaryCondition(is_axially_symmetric, integration_order,
+        new NeumannBoundaryCondition(is_axially_symmetric, integration_order, shapefunction_order,
                                      dof_table, variable_id, component_id,
                                      global_dim, std::move(elements), param));
 }
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
index 0795935b003ab38a154fdc424bd976b9b8bf7019..c956d824a4ae6f12387895172e2b28ec39c54b87 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
@@ -24,7 +24,9 @@ std::unique_ptr<NeumannBoundaryCondition> 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);
 
 }  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
index eab6bc696704388227a495355aa1148046a88ce0..6a3b9e40d168433dc165b1200b885e41a4f6ccbb 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
@@ -17,7 +17,9 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     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)
 {
     DBUG("Constructing RobinBcConfig from config.");
@@ -33,7 +35,7 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     auto const& u_0 = findParameter<double>(u_0_name, parameters, 1);
 
     return std::unique_ptr<RobinBoundaryCondition>(new RobinBoundaryCondition(
-        is_axially_symmetric, integration_order, dof_table, variable_id,
+        is_axially_symmetric, integration_order, shapefunction_order, dof_table, variable_id,
         component_id, global_dim, std::move(elements),
         RobinBoundaryConditionData{alpha, u_0}));
 }
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
index 5c9357f5df290709605dae339e0be8f90d144e3c..64d8f760bd0d4cfe0561aaee31f03abb525ad1c6 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
@@ -40,7 +40,9 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     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);
 
 }  // ProcessLib
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 708fdbd9b8c1bb65e6d698fe9dc2f63ec73538b6..2e7490083dfb1c2c3603d348c1460ef6bb5d2cbc 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -145,7 +145,7 @@ ProcessVariable::createBoundaryConditions(
 
     for (auto& config : _bc_configs)
         bcs.emplace_back(_bc_builder->createBoundaryCondition(
-            config, dof_table, _mesh, variable_id, integration_order, parameters));
+            config, dof_table, _mesh, variable_id, integration_order, _shapefunction_order, parameters));
 
     return bcs;
 }