diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index fec3c23524095b6633d6bc1d873b28a8e01ac8b7..634ffe04fe88817042f028319ea1d22e8104d98c 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -12,6 +12,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshRandomizer.hpp>
 #include <scheme/AcousticSolver.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
@@ -203,6 +204,18 @@ SchemeModule::SchemeModule()
 
                             ));
 
+  this
+    ->_addBuiltinFunction("axis",
+                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
+                            std::shared_ptr<const IBoundaryDescriptor>)>>(
+
+                            [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+                              -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                              return std::make_shared<AxisBoundaryConditionDescriptor>(boundary);
+                            }
+
+                            ));
+
   this
     ->_addBuiltinFunction("symmetry",
                           std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index c6f126b5be3f64451c895c96baae16bf1b26f56d..0743d2d95dda05ca0cb990ade1c8d088c9ed2112 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -2,6 +2,7 @@
 #define MESH_NODE_BOUNDARY_HPP
 
 #include <Kokkos_Vector.hpp>
+#include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
 #include <mesh/ConnectivityMatrix.hpp>
 #include <mesh/IConnectivity.hpp>
@@ -18,6 +19,9 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   Array<const NodeId> m_node_list;
   std::string m_boundary_name;
 
+  template <typename MeshType>
+  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const MeshType& mesh) const;
+
  public:
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
   MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
@@ -76,133 +80,96 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
-    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
-
+  MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
+    : m_boundary_name(ref_edge_list.refId().tagName())
   {
     static_assert(Dimension == MeshType::Dimension);
-  }
-
-  MeshNodeBoundary()          = default;
-  virtual ~MeshNodeBoundary() = default;
-
- protected:
-  MeshNodeBoundary(const MeshNodeBoundary&) = default;
-  MeshNodeBoundary(MeshNodeBoundary&&)      = default;
-};
-
-template <size_t Dimension>
-class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
-{
- public:
-  using Rd = TinyVector<Dimension, double>;
-
- private:
-  const Rd m_outgoing_normal;
-
-  template <typename MeshType>
-  PUGS_INLINE Rd _getNormal(const MeshType& mesh);
+    const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
+    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
 
-  template <typename MeshType>
-  PUGS_INLINE void
-  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
-                       const TinyVector<Dimension, double>& origin,
-                       const double length,
-                       const MeshType& mesh) const
-  {
-    static_assert(MeshType::Dimension == Dimension);
+    auto edge_is_owned = mesh.connectivity().edgeIsOwned();
 
-    const NodeValue<const Rd>& xr = mesh.xr();
+    const Array<const EdgeId>& edge_list = ref_edge_list.list();
 
     bool is_bad = false;
-
-    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
-      const Rd& x = xr[this->m_node_list[r]];
-      if (dot(x - origin, normal) > 1E-13 * length) {
-        is_bad = true;
+    parallel_for(edge_list.size(), [=, &is_bad](int l) {
+      const EdgeId edge_id = edge_list[l];
+      if (edge_is_owned[edge_id]) {
+        const auto& edge_faces             = edge_to_face_matrix[edge_id];
+        bool is_connected_to_boundary_face = false;
+        for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) {
+          const FaceId edge_face_id = edge_faces[i_edge_face];
+          if (face_to_cell_matrix[edge_face_id].size() == 1) {
+            is_connected_to_boundary_face = true;
+          }
+        }
+        if (not is_connected_to_boundary_face) {
+          is_bad = true;
+        }
       }
     });
 
     if (parallel::allReduceOr(is_bad)) {
       std::ostringstream ost;
       ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-          << ": boundary is not flat!";
+          << ": inner edges cannot be used to define mesh boundaries";
       throw NormalError(ost.str());
     }
-  }
 
-  template <typename MeshType>
-  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
+    Kokkos::vector<unsigned int> node_ids;
+    node_ids.reserve(2 * edge_list.size());
+    const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
 
- public:
-  const Rd&
-  outgoingNormal() const
-  {
-    return m_outgoing_normal;
-  }
+    for (size_t l = 0; l < edge_list.size(); ++l) {
+      const EdgeId edge_number = edge_list[l];
+      const auto& edge_nodes   = edge_to_node_matrix[edge_number];
 
-  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+      for (size_t r = 0; r < edge_nodes.size(); ++r) {
+        node_ids.push_back(edge_nodes[r]);
+      }
+    }
+    std::sort(node_ids.begin(), node_ids.end());
+    auto last = std::unique(node_ids.begin(), node_ids.end());
+    node_ids.resize(std::distance(node_ids.begin(), last));
 
-  template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
-  {
-    ;
+    Array<NodeId> node_list(node_ids.size());
+    parallel_for(
+      node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
+    m_node_list = node_list;
   }
 
   template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
+    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
+
   {
-    ;
+    static_assert(Dimension == MeshType::Dimension);
   }
 
-  MeshFlatNodeBoundary()                            = default;
-  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
-  virtual ~MeshFlatNodeBoundary()                   = default;
-};
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 1);
-  using R = TinyVector<1, double>;
-
-  const size_t number_of_bc_nodes = [&]() {
-    size_t number_of_bc_nodes = 0;
-    auto node_is_owned        = mesh.connectivity().nodeIsOwned();
-    for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
-      number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
-    }
-    return parallel::allReduceMax(number_of_bc_nodes);
-  }();
-
-  if (number_of_bc_nodes != 1) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
-        << ": node boundaries in 1D require to have exactly 1 node";
-    throw NormalError(ost.str());
-  }
+  MeshNodeBoundary(const MeshNodeBoundary&) = default;
+  MeshNodeBoundary(MeshNodeBoundary&&)      = default;
 
-  return R{1};
-}
+  MeshNodeBoundary()          = default;
+  virtual ~MeshNodeBoundary() = default;
+};
 
 template <>
 template <typename MeshType>
-PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
+std::array<TinyVector<2>, 2>
+MeshNodeBoundary<2>::_getBounds(const MeshType& mesh) const
 {
   static_assert(MeshType::Dimension == 2);
+
   using R2 = TinyVector<2, double>;
 
   const NodeValue<const R2>& xr = mesh.xr();
 
-  R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
-  R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());
+  std::array<R2, 2> bounds;
+  R2& xmin = bounds[0];
+  R2& xmax = bounds[1];
+
+  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
 
   auto update_xmin = [](const R2& x, R2& xmin) {
     if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
@@ -233,29 +200,16 @@ MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
     }
   }
 
-  if (xmin == xmax) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-        << ": unable to compute normal";
-    throw NormalError(ost.str());
-  }
-
-  R2 dx = xmax - xmin;
-  dx *= 1. / l2Norm(dx);
-
-  const R2 normal(-dx[1], dx[0]);
-
-  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
-
-  return normal;
+  return bounds;
 }
 
 template <>
 template <typename MeshType>
-PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
+std::array<TinyVector<3>, 6>
+MeshNodeBoundary<3>::_getBounds(const MeshType& mesh) const
 {
   static_assert(MeshType::Dimension == 3);
+
   using R3 = TinyVector<3, double>;
 
   auto update_xmin = [](const R3& x, R3& xmin) {
@@ -306,13 +260,21 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
     }
   };
 
-  R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
-  R3 ymin = xmin;
-  R3 zmin = xmin;
+  std::array<R3, 6> bounds;
+  R3& xmin = bounds[0];
+  R3& ymin = bounds[1];
+  R3& zmin = bounds[2];
+  R3& xmax = bounds[3];
+  R3& ymax = bounds[4];
+  R3& zmax = bounds[5];
 
-  R3 xmax = -xmin;
-  R3 ymax = xmax;
-  R3 zmax = xmax;
+  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymin = xmin;
+  zmin = xmin;
+
+  xmax = -xmin;
+  ymax = xmax;
+  zmax = xmax;
 
   const NodeValue<const R3>& xr = mesh.xr();
 
@@ -354,6 +316,157 @@ MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
     }
   }
 
+  return bounds;
+}
+
+template <size_t Dimension>
+class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Rd m_outgoing_normal;
+
+  template <typename MeshType>
+  PUGS_INLINE Rd _getNormal(const MeshType& mesh);
+
+  template <typename MeshType>
+  PUGS_INLINE void
+  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+                       const TinyVector<Dimension, double>& origin,
+                       const double length,
+                       const MeshType& mesh) const
+  {
+    static_assert(MeshType::Dimension == Dimension);
+
+    const NodeValue<const Rd>& xr = mesh.xr();
+
+    bool is_bad = false;
+
+    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+      const Rd& x = xr[this->m_node_list[r]];
+      if (dot(x - origin, normal) > 1E-13 * length) {
+        is_bad = true;
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": boundary is not flat!";
+      throw NormalError(ost.str());
+    }
+  }
+
+  template <typename MeshType>
+  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_outgoing_normal;
+  }
+
+  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+
+  template <typename MeshType>
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+  {
+    ;
+  }
+
+  template <typename MeshType>
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+  {
+    ;
+  }
+
+  MeshFlatNodeBoundary()                            = default;
+  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
+  virtual ~MeshFlatNodeBoundary()                   = default;
+};
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<1, double>
+MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::Dimension == 1);
+  using R = TinyVector<1, double>;
+
+  const size_t number_of_bc_nodes = [&]() {
+    size_t number_of_bc_nodes = 0;
+    auto node_is_owned        = mesh.connectivity().nodeIsOwned();
+    for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
+      number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
+    }
+    return parallel::allReduceMax(number_of_bc_nodes);
+  }();
+
+  if (number_of_bc_nodes != 1) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
+        << ": node boundaries in 1D require to have exactly 1 node";
+    throw NormalError(ost.str());
+  }
+
+  return R{1};
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<2, double>
+MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::Dimension == 2);
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  if (xmin == xmax) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
+  }
+
+  R2 dx = xmax - xmin;
+  dx *= 1. / l2Norm(dx);
+
+  const R2 normal(-dx[1], dx[0]);
+
+  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
+
+  return normal;
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::Dimension == 3);
+
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
+
   const R3 u = xmax - xmin;
   const R3 v = ymax - ymin;
   const R3 w = zmax - zmin;
@@ -530,4 +643,159 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
   }
 }
 
+template <size_t Dimension>
+class MeshLineNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Rd m_direction;
+
+  template <typename MeshType>
+  PUGS_INLINE TinyVector<Dimension> _getDirection(const MeshType&);
+
+  template <typename MeshType>
+  PUGS_INLINE void
+  _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+                       const TinyVector<Dimension, double>& origin,
+                       const double length,
+                       const MeshType& mesh) const
+  {
+    static_assert(MeshType::Dimension == Dimension);
+    using Rdxd = TinyMatrix<Dimension>;
+
+    const NodeValue<const Rd>& xr = mesh.xr();
+
+    Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
+
+    bool is_bad = false;
+    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+      const Rd& x    = xr[this->m_node_list[r]];
+      const Rd delta = x - origin;
+      if (dot(P * delta, direction) > 1E-13 * length) {
+        is_bad = true;
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": boundary is not a line!";
+      throw NormalError(ost.str());
+    }
+  }
+
+ public:
+  const Rd&
+  direction() const
+  {
+    return m_direction;
+  }
+
+  MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  MeshLineNodeBoundary()                            = default;
+  MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary(MeshLineNodeBoundary&&)      = default;
+  virtual ~MeshLineNodeBoundary()                   = default;
+};
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<1>
+MeshLineNodeBoundary<1>::_getDirection(const MeshType&)
+{
+  throw UnexpectedError("MeshLineNodeBoundary::_getDirection make no sense in dimension 1");
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<2>
+MeshLineNodeBoundary<2>::_getDirection(const MeshType& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  if (xmin == xmax) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute direction";
+    throw NormalError(ost.str());
+  }
+
+  R2 direction        = xmax - xmin;
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<3>
+MeshLineNodeBoundary<3>::_getDirection(const MeshType& mesh)
+{
+  static_assert(MeshType::Dimension == 3);
+
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
+
+  const R3 u = xmax - xmin;
+  const R3 v = ymax - ymin;
+  const R3 w = zmax - zmin;
+
+  R3 direction = u;
+  if (dot(v, v) > dot(direction, direction)) {
+    direction = v;
+  }
+
+  if (dot(w, w) > dot(direction, direction)) {
+    direction = w;
+  }
+
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index 36f38e779620cf7844303756a3a54b081da61c9a..1daaff6a04dfbfd4856c8f2d0875bb4cc4549ca5 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -8,6 +8,7 @@
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -27,10 +28,11 @@ class MeshRandomizer
 
   const MeshType& m_given_mesh;
 
-  class SymmetryBoundaryCondition;
+  class AxisBoundaryCondition;
   class FixedBoundaryCondition;
+  class SymmetryBoundaryCondition;
 
-  using BoundaryCondition = std::variant<FixedBoundaryCondition, SymmetryBoundaryCondition>;
+  using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
   BoundaryConditionList m_boundary_condition_list;
@@ -41,6 +43,16 @@ class MeshRandomizer
   {
     BoundaryConditionList bc_list;
 
+    constexpr ItemType EdgeType = [] {
+      if constexpr (Dimension == 3) {
+        return ItemType::edge;
+      } else if constexpr (Dimension == 2) {
+        return ItemType::face;
+      } else {
+        return ItemType::node;
+      }
+    }();
+
     constexpr ItemType FaceType = [] {
       if constexpr (Dimension > 1) {
         return ItemType::face;
@@ -53,6 +65,24 @@ class MeshRandomizer
       bool is_valid_boundary_condition = true;
 
       switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::axis: {
+        const AxisBoundaryConditionDescriptor& axis_bc_descriptor =
+          dynamic_cast<const AxisBoundaryConditionDescriptor&>(*bc_descriptor);
+        for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<EdgeType>();
+             ++i_ref_edge_list) {
+          const auto& ref_edge_list = mesh.connectivity().template refItemList<EdgeType>(i_ref_edge_list);
+          const RefId& ref          = ref_edge_list.refId();
+          if (ref == axis_bc_descriptor.boundaryDescriptor()) {
+            if constexpr (Dimension == 1) {
+              bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_edge_list}});
+            } else {
+              bc_list.emplace_back(AxisBoundaryCondition{MeshLineNodeBoundary<Dimension>(mesh, ref_edge_list)});
+            }
+          }
+        }
+        is_valid_boundary_condition = true;
+        break;
+      }
       case IBoundaryConditionDescriptor::Type::symmetry: {
         const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
           dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -61,8 +91,7 @@ class MeshRandomizer
           const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
           const RefId& ref          = ref_face_list.refId();
           if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-            SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
-            bc_list.push_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
+            bc_list.emplace_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
           }
         }
         is_valid_boundary_condition = true;
@@ -76,11 +105,7 @@ class MeshRandomizer
           const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
           const RefId& ref          = ref_face_list.refId();
           if (ref == fixed_bc_descriptor.boundaryDescriptor()) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
-
-            const auto& node_list = mesh_node_boundary.nodeList();
-
-            bc_list.push_back(FixedBoundaryCondition{node_list});
+            bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_face_list}});
           }
         }
         break;
@@ -121,6 +146,19 @@ class MeshRandomizer
                 shift[node_id] = P * shift[node_id];
               });
 
+          } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+            const Rd& t = bc.direction();
+
+            const Rdxd txt = tensorProduct(t, t);
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                shift[node_id] = txt * shift[node_id];
+              });
+
           } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
             const Array<const NodeId>& node_list = bc.nodeList();
             parallel_for(
@@ -283,52 +321,83 @@ class MeshRandomizer
 };
 
 template <size_t Dimension>
-class MeshRandomizer<Dimension>::SymmetryBoundaryCondition
+class MeshRandomizer<Dimension>::AxisBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
 
  public:
   const Rd&
-  outgoingNormal() const
+  direction() const
   {
-    return m_mesh_flat_node_boundary.outgoingNormal();
+    return m_mesh_line_node_boundary.direction();
   }
 
   const Array<const NodeId>&
   nodeList() const
   {
-    return m_mesh_flat_node_boundary.nodeList();
+    return m_mesh_line_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
-    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+    : m_mesh_line_node_boundary(mesh_line_node_boundary)
   {
     ;
   }
 
-  ~SymmetryBoundaryCondition() = default;
+  ~AxisBoundaryCondition() = default;
 };
 
 template <size_t Dimension>
 class MeshRandomizer<Dimension>::FixedBoundaryCondition
 {
  private:
-  const Array<const NodeId> m_node_list;
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
   nodeList() const
   {
-    return m_node_list;
+    return m_mesh_node_boundary.nodeList();
   }
 
-  FixedBoundaryCondition(const Array<const NodeId>& node_list) : m_node_list{node_list} {}
+  FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
 
   ~FixedBoundaryCondition() = default;
 };
 
+template <size_t Dimension>
+class MeshRandomizer<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_mesh_flat_node_boundary.outgoingNormal();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
 #endif   // MESH_RANDOMIZER_HPP
diff --git a/src/scheme/AxisBoundaryConditionDescriptor.hpp b/src/scheme/AxisBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a03a4bcdd5c2bfbcb23824d3da4e693d4bdbc345
--- /dev/null
+++ b/src/scheme/AxisBoundaryConditionDescriptor.hpp
@@ -0,0 +1,46 @@
+#ifndef AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class AxisBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "axis(" << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::axis;
+  }
+
+  AxisBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+    : m_boundary_descriptor(boundary_descriptor)
+  {
+    ;
+  }
+
+  AxisBoundaryConditionDescriptor(const AxisBoundaryConditionDescriptor&) = delete;
+  AxisBoundaryConditionDescriptor(AxisBoundaryConditionDescriptor&&)      = delete;
+
+  ~AxisBoundaryConditionDescriptor() = default;
+};
+
+#endif   // AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 96989954cc234d5426832b4f8a60a4b3f71f0566..bf92cd450f000db208dfc87c340f2520b9671757 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,6 +8,7 @@ class IBoundaryConditionDescriptor
  public:
   enum class Type
   {
+    axis,
     dirichlet,
     fourier,
     fixed,