diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 311377ea6df090f0d537575ac2b1eba796fea326..fbe44992065514d81e342c683ce80d6d9c1ba088 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -85,15 +85,17 @@ std::unique_ptr<BoundaryCondition> BoundaryConditionBuilder::createBoundaryCondi
         return createNeumannBoundaryCondition(
             config.config,
             getClonedElements(boundary_element_searcher, config.geometry),
-            dof_table, variable_id, config.component_id, integration_order,
-            mesh.getDimension(), parameters);
+            dof_table, variable_id, config.component_id,
+            mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+            parameters);
     }
     else if (type == "Robin") {
         return createRobinBoundaryCondition(
             config.config,
             getClonedElements(boundary_element_searcher, config.geometry),
-            dof_table, variable_id, config.component_id, integration_order,
-            mesh.getDimension(), parameters);
+            dof_table, variable_id, config.component_id,
+            mesh.isAxiallySymmetric(), integration_order, mesh.getDimension(),
+            parameters);
     }
     else
     {
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index cc81024ed1fe729d59bd6e7b3e08d48a63b474fb..515b653b8592c1332358f0e91c71e93cc1e2a89b 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -24,11 +24,12 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
         typename std::enable_if<
             std::is_same<typename std::decay<BoundaryConditionData>::type,
                          typename std::decay<Data>::type>::value,
-            unsigned const>::type integration_order,
+            bool>::type is_axially_symmetric,
+        unsigned const integration_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
-        unsigned const global_dim,
-        std::vector<MeshLib::Element*>&& elements, Data&& data)
+        unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+        Data&& data)
     : _data(std::forward<Data>(data)),
       _elements(std::move(elements)),
       _integration_order(integration_order)
@@ -55,7 +56,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
 
     createLocalAssemblers<LocalAssemblerImplementation>(
         global_dim, _elements, *_dof_table_boundary, _local_assemblers,
-        false /* TODO */, _integration_order, _data);
+        is_axially_symmetric, _integration_order, _data);
 }
 
 template <typename BoundaryConditionData,
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index f40fc809009a51d35509888a8ce32c9948623017..5edcf191a7302fc8a42439aa265a22fa8700a902 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -31,11 +31,12 @@ public:
         typename std::enable_if<
             std::is_same<typename std::decay<BoundaryConditionData>::type,
                          typename std::decay<Data>::type>::value,
-            unsigned const>::type integration_order,
+            bool>::type is_axially_symmetric,
+        unsigned const integration_order,
         NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
         int const variable_id, int const component_id,
-        unsigned const global_dim,
-        std::vector<MeshLib::Element*>&& elements, Data&& data);
+        unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+        Data&& data);
 
     ~GenericNaturalBoundaryCondition() override;
 
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
index 448ef07b5354b2ae2a07c555f30b54d8a20ef8f5..79601405d89f4dc9446a89cdb636d37a45a83a74 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
@@ -12,13 +12,12 @@
 
 namespace ProcessLib
 {
-std::unique_ptr<NeumannBoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     BaseLib::ConfigTree const& config,
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, unsigned const integration_order,
-    unsigned const global_dim,
+    int const component_id, bool is_axially_symmetric,
+    unsigned const integration_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Constructing Neumann BC from config.");
@@ -32,9 +31,9 @@ createNeumannBoundaryCondition(
     auto const& param = findParameter<double>(param_name, parameters, 1);
 
     return std::unique_ptr<NeumannBoundaryCondition>(
-        new NeumannBoundaryCondition(
-            integration_order, dof_table, variable_id, component_id,
-            global_dim, std::move(elements), param));
+        new NeumannBoundaryCondition(is_axially_symmetric, integration_order,
+                                     dof_table, variable_id, component_id,
+                                     global_dim, std::move(elements), param));
 }
 
 }  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
index 53713fb255bda0488943968e0978deae81afa799..0795935b003ab38a154fdc424bd976b9b8bf7019 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
@@ -19,13 +19,12 @@ namespace ProcessLib
 using NeumannBoundaryCondition = GenericNaturalBoundaryCondition<
     Parameter<double> const&, NeumannBoundaryConditionLocalAssembler>;
 
-std::unique_ptr<NeumannBoundaryCondition>
-createNeumannBoundaryCondition(
+std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     BaseLib::ConfigTree const& config,
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, unsigned const integration_order,
-    unsigned const global_dim,
+    int const component_id, bool is_axially_symmetric,
+    unsigned const integration_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 df8c74183ae35cc60463ae80d0db574643ba5ef1..eab6bc696704388227a495355aa1148046a88ce0 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
@@ -16,8 +16,8 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     BaseLib::ConfigTree const& config,
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, unsigned const integration_order,
-    unsigned const global_dim,
+    int const component_id, bool is_axially_symmetric,
+    unsigned const integration_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters)
 {
     DBUG("Constructing RobinBcConfig from config.");
@@ -32,11 +32,10 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     auto const& alpha = findParameter<double>(alpha_name, parameters, 1);
     auto const& u_0 = findParameter<double>(u_0_name, parameters, 1);
 
-    return std::unique_ptr<RobinBoundaryCondition>(
-        new RobinBoundaryCondition(
-            integration_order, dof_table, variable_id, component_id, global_dim,
-            std::move(elements),
-            RobinBoundaryConditionData{alpha, u_0}));
+    return std::unique_ptr<RobinBoundaryCondition>(new RobinBoundaryCondition(
+        is_axially_symmetric, integration_order, dof_table, variable_id,
+        component_id, global_dim, std::move(elements),
+        RobinBoundaryConditionData{alpha, u_0}));
 }
 
 }  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
index 2b67f531c018b70a476c38bc8fe26137dea3431c..5c9357f5df290709605dae339e0be8f90d144e3c 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.h
@@ -39,8 +39,8 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     BaseLib::ConfigTree const& config,
     std::vector<MeshLib::Element*>&& elements,
     NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id, unsigned const integration_order,
-    unsigned const global_dim,
+    int const component_id, bool is_axially_symmetric,
+    unsigned const integration_order, unsigned const global_dim,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
 }  // ProcessLib