diff --git a/include/bout/boundary_common.hxx b/include/bout/boundary_common.hxx
new file mode 100644
index 0000000000..6d08f2c9c0
--- /dev/null
+++ b/include/bout/boundary_common.hxx
@@ -0,0 +1,57 @@
+#pragma once
+
+#include <bout/bout_types.hxx>
+
+#include <cstdint>
+
+namespace bout {
+namespace boundary {
+/// Types of free boundary condition
+/// ================== ===================================================
+/// Name               Description
+/// ================== ===================================================
+/// ``limited``        use exponential if decreasing, otherwise Neumanm
+/// ``exponential``    use exponential extrapolation
+/// ``linear``         use linear extrapolation
+/// ================== ===================================================
+enum class BoundaryFreeExtrapolation : std::uint8_t { limited, exponential, linear };
+
+// Limited free gradient of log of a quantity
+// This ensures that the guard cell values remain positive
+// while also ensuring that the quantity never increases
+//
+//  fm  fc | fp
+//         ^ boundary
+//
+// exp( 2*log(fc) - log(fm) )
+inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc,
+                               bout::boundary::BoundaryFreeExtrapolation mode) {
+  if ((fm < fc) && (mode == bout::boundary::BoundaryFreeExtrapolation::limited)) {
+    return fc; // Neumann rather than increasing into boundary
+  }
+  if (fm < 1e-10) {
+    return fc; // Low / no density condition
+  }
+
+  BoutReal fp = 0;
+  switch (mode) {
+  case bout::boundary::BoundaryFreeExtrapolation::limited:
+  case bout::boundary::BoundaryFreeExtrapolation::exponential:
+    fp = fc * fc / fm; // Exponential
+    break;
+  case bout::boundary::BoundaryFreeExtrapolation::linear:
+    fp = (2.0 * fc) - fm; // Linear
+    break;
+  }
+
+#if CHECKLEVEL >= 2
+  if (!std::isfinite(fp)) {
+    throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp);
+  }
+#endif
+
+  return fp;
+}
+
+} // namespace boundary
+} // namespace bout
diff --git a/include/bout/boundary_factory.hxx b/include/bout/boundary_factory.hxx
index 5f1f6e06a6..17807138be 100644
--- a/include/bout/boundary_factory.hxx
+++ b/include/bout/boundary_factory.hxx
@@ -13,17 +13,18 @@ class BoundaryRegion;
 class BoundaryModifier;
 
 #include <map>
+#include <memory>
 #include <string>
 
 /// Create BoundaryOp objects on demand
 /*!
  * This implements a simple string parser, used to match boundary condition
  * names like "dirichlet" with a BoundaryOp object.
- * 
+ *
  * Modifiers: Simple modifications of boundary conditions can be performed,
  * for example transforming the coordinate system
- * 
- * This is a singleton, so only one instance can exist. This is 
+ *
+ * This is a singleton, so only one instance can exist. This is
  * enforced by making the constructor private, and having a getInstance() method
  * to return a pointer to the only instance.
  *
@@ -33,7 +34,7 @@ class BoundaryModifier;
  * Boundaries are defined as classes which inherit from BoundaryOp
  * These define a clone() function which creates a new BoundaryOp,
  * given a list of arguments. See boundary_standard.hxx for examples.
- * 
+ *
  *     class MyBoundary : public BoundaryOp {
  *      public:
  *       BoundaryOp* clone(BoundaryRegion *region, const list<string> &args) {
@@ -43,23 +44,23 @@ class BoundaryModifier;
  *       void apply(Field2D &f);
  *       void apply(Field3D &f);
  *     };
- * 
+ *
  * The singleton instance of BoundaryFactory from getInstance():
- * 
+ *
  *     BoundaryFactory* bf = BoundaryFactory::getInstance();
  *
  * New boundary types can be added to the BoundaryFactory
  *
  *     bf->add(new MyBoundary, "myboundary");
- * 
+ *
  *
  * Subsequent calls to create() or createFromOptions() can make use
  * of the boundary type "myboundary".
  *
  * BoundaryOpBase *bndry = bf->create("myboundary()", new BoundaryRegionXOut("xout", 0, 10, localmesh));
- * 
+ *
  * where the region is defined in boundary_region.hxx
- * 
+ *
  */
 class BoundaryFactory {
 public:
@@ -70,13 +71,16 @@ public:
   static void cleanup(); ///< Frees all memory
 
   /// Create a boundary operation object
-  BoundaryOpBase* create(const std::string& name, BoundaryRegionBase* region);
-  BoundaryOpBase* create(const char* name, BoundaryRegionBase* region);
+  BoundaryOpBase* create(const std::string& name,
+                         const std::shared_ptr<BoundaryRegionBase>& region);
+  BoundaryOpBase* create(const char* name,
+                         const std::shared_ptr<BoundaryRegionBase>& region);
 
   /// Create a boundary object using the options file
   BoundaryOpBase* createFromOptions(const std::string& varname,
-                                    BoundaryRegionBase* region);
-  BoundaryOpBase* createFromOptions(const char* varname, BoundaryRegionBase* region);
+                                    const std::shared_ptr<BoundaryRegionBase>& region);
+  BoundaryOpBase* createFromOptions(const char* varname,
+                                    const std::shared_ptr<BoundaryRegionBase>& region);
 
   /*!
    * Add available boundary conditions and modifiers
diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx
index bf03e15636..b2d2b9fb99 100644
--- a/include/bout/boundary_iterator.hxx
+++ b/include/bout/boundary_iterator.hxx
@@ -38,122 +38,6 @@ protected:
 public:
   BoundaryRegionIter& operator*() { return *this; }
 
-  void dirichlet_o2(Field3D& f, BoutReal value) const {
-    ynext(f) = bout::parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value);
-  }
-
-  BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); }
-
-  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
-    return (f[ind()] * 3 - yprev(f)) * 0.5;
-  }
-
-  BoutReal extrapolate_next_o2(const Field3D& f) const {
-    return (2 * f[ind()]) - yprev(f);
-  }
-
-  BoutReal
-  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return (2 * f(0, ind())) - f(0, ind().yp(-_by).xp(-_bx));
-  }
-
-  BoutReal interpolate_sheath_o2(const Field3D& f) const {
-    return (f[ind()] + ynext(f)) * 0.5;
-  }
-
-  BoutReal
-  interpolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return (f(0, ind()) + f(0, ind().yp(-_by).xp(-_bx))) * 0.5;
-  }
-
-  BoutReal
-  extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-_by).xp(-_bx)));
-  }
-
-  BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const {
-    const BoutReal fac =
-        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
-    const BoutReal val = ythis(f);
-    const BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac;
-    return 0.5 * (val + next);
-  }
-
-  void set_free(Field3D& f, SheathLimitMode mode) const {
-    const BoutReal fac =
-        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
-    BoutReal val = ythis(f);
-    if (mode == SheathLimitMode::linear_free) {
-      for (int i = 1; i <= localmesh->ystart; ++i) {
-        val += fac;
-        f[ind().yp(_by * i).xp(_bx * i)] = val;
-      }
-    } else {
-      for (int i = 1; i <= localmesh->ystart; ++i) {
-        val *= fac;
-        f[ind().yp(_by * i).xp(_bx * i)] = val;
-      }
-    }
-  }
-
-  void limitFree(Field3D& f) const {
-    const BoutReal fac =
-        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f));
-    BoutReal val = ythis(f);
-    for (int i = 1; i <= localmesh->ystart; ++i) {
-      val *= fac;
-      f[ind().yp(_by * i).xp(_bx * i)] = val;
-    }
-  }
-
-  bool is_lower() const {
-    ASSERT2(_bx == 0);
-    return _by == -1;
-  }
-
-  void neumann_o1(Field3D& f, BoutReal grad) const {
-    BoutReal val = ythis(f);
-    for (int i = 1; i <= localmesh->ystart; ++i) {
-      val += grad;
-      f[ind().yp(_by * i).xp(_bx * i)] = val;
-    }
-  }
-
-  void neumann_o2(Field3D& f, BoutReal grad) const {
-    BoutReal val = yprev(f) + grad;
-    for (int i = 1; i <= localmesh->ystart; ++i) {
-      val += grad;
-      f[ind().yp(_by * i).xp(_bx * i)] = val;
-    }
-  }
-
-  void limit_at_least(Field3D& f, BoutReal value) const {
-    ynext(f) = std::max(ynext(f), value);
-  }
-
-  BoutReal& ynext(Field3D& f) const { return f[ind().yp(_by).xp(_bx)]; }
-  const BoutReal& ynext(const Field3D& f) const { return f[ind().yp(_by).xp(_bx)]; }
-  BoutReal& yprev(Field3D& f) const { return f[ind().yp(-_by).xp(-_bx)]; }
-  const BoutReal& yprev(const Field3D& f) const { return f[ind().yp(-_by).xp(-_bx)]; }
-  BoutReal& ythis(Field3D& f) const { return f[ind()]; }
-  const BoutReal& ythis(const Field3D& f) const { return f[ind()]; }
-
-  void setYPrevIfValid(Field3D& f, BoutReal val) const { yprev(f) = val; }
-  void setAll(Field3D& f, const BoutReal val) const {
-    for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) {
-      f[ind().yp(_by * i).xp(_bx * i)] = val;
-    }
-  }
-
-  static int abs_offset() { return 1; }
-
-  BoutReal& ynext(Field2D& f) const { return f[ind().yp(_by).xp(_bx)]; }
-  const BoutReal& ynext(const Field2D& f) const { return f[ind().yp(_by).xp(_bx)]; }
-  BoutReal& yprev(Field2D& f) const { return f[ind().yp(-_by).xp(-_bx)]; }
-  const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-_by).xp(-_bx)]; }
-
-  int dir() const { return _dir; }
-
 protected:
   int x() const { return _x; }
   int y() const { return _y; }
diff --git a/include/bout/boundary_region_iter.hxx b/include/bout/boundary_region_iter.hxx
index 242306c04f..6b6cd5b6ff 100644
--- a/include/bout/boundary_region_iter.hxx
+++ b/include/bout/boundary_region_iter.hxx
@@ -5,11 +5,13 @@
 #include <cstddef>
 #include <cstdint>
 #include <functional>
+#include <memory>
 #include <string>
 #include <type_traits>
 #include <vector>
 
 #include "bout/assert.hxx"
+#include "bout/boundary_common.hxx"
 #include "bout/bout_types.hxx"
 #include "bout/field_data.hxx"
 #include "bout/utils.hxx"
@@ -34,17 +36,6 @@ enum class BndryType : std::int8_t {
   num
 };
 
-/// Types of free boundary condition
-/// ================== ===================================================
-/// Name               Description
-/// ================== ===================================================
-/// ``limited``        use exponential if decreasing, otherwise Neumanm
-/// ``exponential``    use exponential extrapolation
-/// ``linear``         use linear extrapolation
-/// ================== ===================================================
-enum class BoundaryFreeExtrapolation : std::int8_t { limited, exponential, linear };
-//BOUT_ENUM_CLASS(BoundaryFreeExtrapolation, limited, exponential, linear);
-
 template <typename impl>
 class BoundaryRegionIterBase {
   BoundaryRegionIterBase() = default;
@@ -156,23 +147,20 @@ public:
   /// off = 1 is the last point in the domain
   /// off = 2 is the second to last point in the domain
   template <bool check = true>
-  BoutReal& getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& func,
-                  int off) const {
+  BoutReal getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& func,
+                 int off) const {
     return static_cast<const impl*>(this)->template _getAt<check>(func, off);
   }
   /// Get the first point in the boundary
-  const BoutReal&
-  next(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
+  BoutReal next(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
     return static_cast<const impl*>(this)->_getAt(func, 0);
   }
   /// Get the last point in the domain
-  const BoutReal&
-  current(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
+  BoutReal current(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
     return static_cast<const impl*>(this)->_getAt(func, 1);
   }
   /// Get the second to last point in the domain - this may not be valid and thus throw
-  const BoutReal&
-  prev(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
+  BoutReal prev(const std::function<BoutReal(int yoffset, Ind3D ind)>& func) const {
     return static_cast<const impl*>(this)->_getAt(func, 2);
   }
 
@@ -528,8 +516,8 @@ public:
     return f.ynext(_off)[_ind().yp(_off)];
   }
   template <bool check = true>
-  BoutReal getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
-                 int off) const {
+  BoutReal _getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
+                  int off) const {
     if constexpr (check) {
       ASSERT3(valid() > -off - 2);
     }
@@ -571,7 +559,7 @@ public:
       location = dir == 1 ? BNDRY_XOUT : BNDRY_XIN;
     } else {
       this->isY = true;
-      location = dir == 1 ? BNDRY_YDOWN : BNDRY_YUP;
+      location = dir == 1 ? BNDRY_YUP : BNDRY_YDOWN;
     }
   }
   int dir() { return _dir; }
@@ -648,8 +636,8 @@ public:
     }
   }
   template <bool check = true>
-  BoutReal getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
-                 int off) const {
+  BoutReal _getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
+                  int off) const {
     if constexpr (check) {
       ASSERT3(_valid() > -off - 2);
     }
@@ -694,9 +682,9 @@ using BoundaryRegionY = BoundaryRegionXY<false>;
 using BoundaryRegionIterX = BoundaryRegionIterXY<true>;
 using BoundaryRegionIterY = BoundaryRegionIterXY<false>;
 
-inline BoundaryRegionX* NewBoundaryRegionXIn(const std::string& name, int ymin, int ymax,
-                                             Mesh* mesh) {
-  auto* pointer = new BoundaryRegionX(
+inline std::shared_ptr<BoundaryRegionX>
+NewBoundaryRegionXIn(const std::string& name, int ymin, int ymax, Mesh* mesh) {
+  auto pointer = std::make_shared<BoundaryRegionX>(
       name, -1, mesh,
       Region<Ind3D>(mesh->xstart, mesh->xstart, ymin, ymax, mesh->zstart, mesh->zend,
                     mesh->LocalNy, mesh->LocalNz, mesh->maxregionblocksize));
@@ -704,9 +692,9 @@ inline BoundaryRegionX* NewBoundaryRegionXIn(const std::string& name, int ymin,
   return pointer;
 }
 
-inline BoundaryRegionX* NewBoundaryRegionXOut(const std::string& name, int ymin, int ymax,
-                                              Mesh* mesh) {
-  auto* pointer = new BoundaryRegionX(
+inline std::shared_ptr<BoundaryRegionX>
+NewBoundaryRegionXOut(const std::string& name, int ymin, int ymax, Mesh* mesh) {
+  auto pointer = std::make_shared<BoundaryRegionX>(
       name, 1, mesh,
       Region<Ind3D>(mesh->xend, mesh->xend, ymin, ymax, mesh->zstart, mesh->zend,
                     mesh->LocalNy, mesh->LocalNz, mesh->maxregionblocksize));
@@ -714,20 +702,20 @@ inline BoundaryRegionX* NewBoundaryRegionXOut(const std::string& name, int ymin,
   return pointer;
 }
 
-inline BoundaryRegionY* NewBoundaryRegionYUp(const std::string& name, int xmin, int xmax,
-                                             Mesh* mesh) {
-  auto* pointer = new BoundaryRegionY(
-      name, -1, mesh,
+inline std::shared_ptr<BoundaryRegionY>
+NewBoundaryRegionYUp(const std::string& name, int xmin, int xmax, Mesh* mesh) {
+  auto pointer = std::make_shared<BoundaryRegionY>(
+      name, 1, mesh,
       Region<Ind3D>(xmin, xmax, mesh->yend, mesh->yend, mesh->zstart, mesh->zend,
                     mesh->LocalNy, mesh->LocalNz, mesh->maxregionblocksize));
   pointer->legacy = new ::BoundaryRegionYUp(name, xmin, xmax, mesh);
   return pointer;
 }
 
-inline BoundaryRegionY* NewBoundaryRegionYDown(const std::string& name, int xmin,
-                                               int xmax, Mesh* mesh) {
-  auto* pointer = new BoundaryRegionY(
-      name, 1, mesh,
+inline std::shared_ptr<BoundaryRegionY>
+NewBoundaryRegionYDown(const std::string& name, int xmin, int xmax, Mesh* mesh) {
+  auto pointer = std::make_shared<BoundaryRegionY>(
+      name, -1, mesh,
       Region<Ind3D>(xmin, xmax, mesh->ystart, mesh->ystart, mesh->zstart, mesh->zend,
                     mesh->LocalNy, mesh->LocalNz, mesh->maxregionblocksize));
   pointer->legacy = new ::BoundaryRegionYDown(name, xmin, xmax, mesh);
@@ -735,17 +723,18 @@ inline BoundaryRegionY* NewBoundaryRegionYDown(const std::string& name, int xmin
 }
 
 template <class Func>
-void iter_boundary(const BoundaryRegionBase* bndrybase, const Func& func) {
+void iter_boundary(const std::shared_ptr<const BoundaryRegionBase>& bndrybase,
+                   const Func& func) {
   if (bndrybase->isX) {
-    const auto* const bndry = dynamic_cast<const BoundaryRegionX*>(bndrybase);
+    const auto* const bndry = dynamic_cast<const BoundaryRegionX*>(bndrybase.get());
     return iter_boundary(*bndry, func);
   }
   if (bndrybase->isY) {
-    const auto* const bndry = dynamic_cast<const BoundaryRegionY*>(bndrybase);
+    const auto* const bndry = dynamic_cast<const BoundaryRegionY*>(bndrybase.get());
     return iter_boundary(*bndry, func);
   }
   if (bndrybase->isParallel) {
-    const auto* const bndry = dynamic_cast<const BoundaryRegionFCI*>(bndrybase);
+    const auto* const bndry = dynamic_cast<const BoundaryRegionFCI*>(bndrybase.get());
     return iter_boundary(*bndry, func);
   }
   throw BoutException("{} is of unknown type - probably a legacy iterator",
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 56baea3149..84f9ec36c6 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -486,7 +486,7 @@ public:
   // Boundary regions
 
   /// Return a vector containing all the boundary regions on this processor
-  virtual std::vector<BoundaryRegionBase*> getBoundaries() const = 0;
+  virtual std::vector<std::shared_ptr<BoundaryRegionBase>> getBoundaries() const = 0;
 
   /// Get the set of all possible boundaries in this configuration
   virtual std::set<std::string> getPossibleBoundaries() const {
@@ -494,7 +494,7 @@ public:
   };
 
   /// Add a boundary region to this processor
-  virtual void addBoundary(BoundaryRegionBase* UNUSED(bndry)) {
+  virtual void addBoundary(std::shared_ptr<BoundaryRegionBase> UNUSED(bndry)) {
     throw BoutException("This has never been implemented");
   };
 
diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx
index d6fac3c796..2d17288a3d 100644
--- a/include/bout/parallel_boundary_region.hxx
+++ b/include/bout/parallel_boundary_region.hxx
@@ -26,8 +26,6 @@
  *
  */
 
-BOUT_ENUM_CLASS(SheathLimitMode, limit_free, exponential_free, linear_free);
-
 namespace bout {
 namespace parallel_boundary_region {
 
@@ -65,58 +63,6 @@ using IndicesVec = std::vector<Indices>;
 using IndicesIter = IndicesVec::iterator;
 using IndicesIterConst = IndicesVec::const_iterator;
 
-/// Limited free gradient of log of a quantity
-/// This ensures that the guard cell values remain positive
-/// while also ensuring that the quantity never increases
-///
-///  fm  fc | fp
-///         ^ boundary
-///
-/// exp( 2*log(fc) - log(fm) )
-inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc, SheathLimitMode mode) {
-  if ((fm < fc) && (mode == SheathLimitMode::limit_free)) {
-    return fc; // Neumann rather than increasing into boundary
-  }
-  if (fm < 1e-10) {
-    return fc; // Low / no density condition
-  }
-
-  BoutReal fp = 0;
-  switch (mode) {
-  case SheathLimitMode::limit_free:
-  case SheathLimitMode::exponential_free:
-    fp = SQ(fc) / fm; // Exponential
-    break;
-  case SheathLimitMode::linear_free:
-    fp = (2.0 * fc) - fm; // Linear
-    break;
-  }
-
-#if CHECKLEVEL >= 2
-  if (!std::isfinite(fp)) {
-    throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp);
-  }
-#endif
-
-  return fp;
-}
-
-inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) {
-  if (fm < fc) {
-    return 1; // Neumann rather than increasing into boundary
-  }
-  if (fm < 1e-10) {
-    return 1; // Low / no density condition
-  }
-  BoutReal fp = fc / fm;
-#if CHECKLEVEL >= 2
-  if (!std::isfinite(fp)) {
-    throw BoutException("SheathBoundaryParallel limitFree: {}, {} -> {}", fm, fc, fp);
-  }
-#endif
-  return fp;
-}
-
 template <class IndicesVec, class IndicesIter>
 class BoundaryRegionParIterBase {
 
@@ -144,258 +90,6 @@ public:
   // setter
   void setValid(signed char valid) { bndry_position->valid = valid; }
 
-  // extrapolate a given point to the boundary
-  BoutReal extrapolate_sheath_o1(const Field3D& f) const { return ythis(f); }
-  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return extrapolate_sheath_o1(f);
-    }
-    return ythis(f) * (1 + length(f.getLocation())) - yprev(f) * length(f.getLocation());
-  }
-  BoutReal
-  extrapolate_sheath_o1(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return ythis(f);
-  }
-  BoutReal extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
-                                 CELL_LOC loc = CELL_CENTRE) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return extrapolate_sheath_o1(f);
-    }
-    return ythis(f) * (1 + length(loc)) - yprev(f) * length(loc);
-  }
-
-  BoutReal interpolate_sheath_o2(const Field3D& f) const {
-    return ythis(f) * (1 - length(f.getLocation())) + ynext(f) * length(f.getLocation());
-  }
-  BoutReal interpolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
-                                 CELL_LOC loc = CELL_CENTRE) const {
-    return ythis(f) * (1 - length(loc)) + ynext(f) * length(loc);
-  }
-
-  BoutReal extrapolate_next_o1(const Field3D& f) const { return ythis(f); }
-  BoutReal extrapolate_next_o2(const Field3D& f) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return extrapolate_next_o1(f);
-    }
-    return ythis(f) * 2 - yprev(f);
-  }
-
-  BoutReal
-  extrapolate_next_o1(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return ythis(f);
-  }
-  BoutReal
-  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return extrapolate_next_o1(f);
-    }
-    return ythis(f) * 2 - yprev(f);
-  }
-
-  // extrapolate the gradient into the boundary
-  BoutReal extrapolate_grad_o1([[maybe_unused]] const Field3D& f) const { return 0; }
-  BoutReal extrapolate_grad_o2(const Field3D& f) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return extrapolate_grad_o1(f);
-    }
-    return ythis(f) - ynext(f);
-  }
-
-  BoundaryRegionParIterBase& operator*() { return *this; }
-
-  BoundaryRegionParIterBase& operator++() {
-    ++bndry_position;
-    return *this;
-  }
-
-  bool operator!=(const BoundaryRegionParIterBase& rhs) {
-    return bndry_position != rhs.bndry_position;
-  }
-
-#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i)
-  // dirichlet boundary code
-  void dirichlet_o1(Field3D& f, BoutReal value) const {
-    ITER() { getAt(f, i) = value; }
-  }
-
-  void dirichlet_o2(Field3D& f, BoutReal value) const {
-    if (length(f.getLocation()) < small_value) {
-      return dirichlet_o1(f, value);
-    }
-    ITER() {
-      getAt(f, i) = parallel_stencil::dirichlet_o2(
-          i + 1, ythis(f), i + 1 - length(f.getLocation()), value);
-    }
-  }
-
-  void dirichlet_o3(Field3D& f, BoutReal value) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return dirichlet_o2(f, value);
-    }
-    if (length(f.getLocation()) < small_value) {
-      ITER() {
-        getAt(f, i) = parallel_stencil::dirichlet_o2(
-            i + 2, yprev(f), i + 1 - length(f.getLocation()), value);
-      }
-    } else {
-      ITER() {
-        getAt(f, i) = parallel_stencil::dirichlet_o3(
-            i + 2, yprev(f), i + 1, ythis(f), i + 1 - length(f.getLocation()), value);
-      }
-    }
-  }
-
-  void limit_at_least(Field3D& f, BoutReal value) const {
-    ITER() {
-      if (getAt(f, i) < value) {
-        getAt(f, i) = value;
-      }
-    }
-  }
-
-  bool is_lower() const { return _dir == -1; }
-
-  // NB: value needs to be scaled by dy
-  // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil.
-  // But in general we do not, and thus for normal C2 stencils, this is 1st order.
-  void neumann_o1(Field3D& f, BoutReal value) const {
-    ITER() { getAt(f, i) = ythis(f) + value * (i + 1); }
-  }
-
-  // NB: value needs to be scaled by dy
-  void neumann_o2(Field3D& f, BoutReal value) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return neumann_o1(f, value);
-    }
-    ITER() { getAt(f, i) = yprev(f) + (2 + i) * value; }
-  }
-
-  // NB: value needs to be scaled by dy
-  void neumann_o3(Field3D& f, BoutReal value) const {
-    ASSERT3(valid() >= 0);
-    if (valid() < 1) {
-      return neumann_o2(f, value);
-    }
-    ITER() {
-      getAt(f, i) = parallel_stencil::neumann_o3(i + 1 - length(f.getLocation()), value,
-                                                 i + 1, ythis(f), 2, yprev(f));
-    }
-  }
-
-  // extrapolate into the boundary using only monotonic decreasing values.
-  // f needs to be positive
-  void limitFree(Field3D& f) const {
-    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f)) : 1;
-    auto val = ythis(f);
-    ITER() {
-      val *= fac;
-      getAt(f, i) = val;
-    }
-  }
-
-  BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const {
-    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode)
-                                 : (mode == SheathLimitMode::linear_free ? 0 : 1);
-    auto val = ythis(f);
-    BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac;
-    return val * length(f.getLocation()) + next * (1 - length(f.getLocation()));
-  }
-
-  void set_free(Field3D& f, SheathLimitMode mode) const {
-    const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode)
-                                 : (mode == SheathLimitMode::linear_free ? 0 : 1);
-    auto val = ythis(f);
-    if (mode == SheathLimitMode::linear_free) {
-      ITER() {
-        val += fac;
-        getAt(f, i) = val;
-      }
-    } else {
-      ITER() {
-        val *= fac;
-        getAt(f, i) = val;
-      }
-    }
-  }
-
-  void setAll(Field3D& f, const BoutReal val) const {
-    for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) {
-      f.ynext(i)[ind().yp(i)] = val;
-    }
-  }
-
-  template <bool check = true>
-  BoutReal& getAt(Field3D& f, int off) const {
-    ASSERT3(f.hasParallelSlices());
-    if constexpr (check) {
-      ASSERT3(valid() > -off - 2);
-    }
-    auto _off = offset() + off * _dir;
-    return f.ynext(_off)[ind().yp(_off)];
-  }
-  template <bool check = true>
-  const BoutReal& getAt(const Field3D& f, int off) const {
-    ASSERT3(f.hasParallelSlices());
-    if constexpr (check) {
-      ASSERT3(valid() > -off - 2);
-    }
-    auto _off = offset() + off * _dir;
-    return f.ynext(_off)[ind().yp(_off)];
-  }
-
-  const BoutReal& ynext(const Field3D& f) const { return getAt(f, 0); }
-  BoutReal& ynext(Field3D& f) const { return getAt(f, 0); }
-  const BoutReal& ythis(const Field3D& f) const { return getAt(f, -1); }
-  BoutReal& ythis(Field3D& f) const { return getAt(f, -1); }
-  const BoutReal& yprev(const Field3D& f) const { return getAt(f, -2); }
-  BoutReal& yprev(Field3D& f) const { return getAt(f, -2); }
-
-  template <bool check = true>
-  BoutReal getAt(const std::function<BoutReal(int yoffset, Ind3D ind)>& f,
-                 int off) const {
-    if constexpr (check) {
-      ASSERT3(valid() > -off - 2);
-    }
-    auto _off = offset() + off * _dir;
-    return f(_off, ind().yp(_off));
-  }
-  BoutReal ynext(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return getAt(f, 0);
-  }
-  BoutReal ythis(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return getAt(f, -1);
-  }
-  BoutReal yprev(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
-    return getAt(f, -2);
-  }
-
-  void setYPrevIfValid(Field3D& f, BoutReal val) const {
-    if (valid() > 0) {
-      yprev(f) = val;
-    }
-  }
-
-#if not(BOUT_USE_METRIC_3D)
-  const BoutReal& ynext(const Field2D& f) const { return f.ynext(_dir)[ind().yp(_dir)]; }
-  BoutReal& ynext(Field2D& f) const { return f.ynext(_dir)[ind().yp(_dir)]; }
-
-  const BoutReal& yprev(const Field2D& f) const {
-    ASSERT3(valid() > 0);
-    return f.ynext(-_dir)[ind().yp(-_dir)];
-  }
-  BoutReal& yprev(Field2D& f) const {
-    ASSERT3(valid() > 0);
-    return f.ynext(-_dir)[ind().yp(-_dir)];
-  }
-#endif
-
 private:
   const IndicesVec& bndry_points;
   IndicesIter bndry_position;
diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx
index b1aec6bc2e..bb9f52bbda 100644
--- a/include/bout/yboundary_regions.hxx
+++ b/include/bout/yboundary_regions.hxx
@@ -1,9 +1,9 @@
 #pragma once
 
-#include <bout/boundary_region_iter.hxx>
+#include "bout/boundary_region_iter.hxx"
 
-#include "./boundary_iterator.hxx"
 #include "bout/assert.hxx"
+#include "bout/boundary_iterator.hxx"
 #include "bout/boutexception.hxx"
 #include "bout/field_data.hxx"
 #include "bout/globals.hxx"
@@ -24,17 +24,12 @@
 
 class YBoundary {
 public:
+  /// Iterate over the boundary.
+  /// This function takes a lamda / templated function, that applies the boundary on the given point.
+  /// The function must take a `auto& point` as argument.
+  /// See also the documentation at ../../manual/sphinx/user_docs/boundary_options.rst
   template <class Func>
-  void iter_regions(const Func& func) {
-    for (auto& region : boundary_regions) {
-      func(*region);
-    }
-    for (auto& region : boundary_regions_par) {
-      func(*region);
-    }
-  }
-  template <class Func>
-  void iter_points(const Func& func) {
+  void iter(const Func& func) {
     iter_regions([&](auto& region) {
       for (auto& point : region) {
         func(point);
@@ -42,11 +37,6 @@ public:
     });
   }
 
-  template <class Func>
-  void iter(const Func& func) {
-    return iter_points(func);
-  }
-
   YBoundary(YBndryType type, Options* options_ptr, const Mesh& mesh) {
     bool lower_y = true;
     bool upper_y = true;
@@ -94,19 +84,18 @@ public:
         }
       }
     } else {
-      if (lower_y) {
-        boundary_regions.push_back(
-            std::make_shared<NewBoundaryRegionY>(mesh, true, mesh.iterateBndryLowerY()));
-      }
-      if (upper_y) {
-        boundary_regions.push_back(
-            std::make_shared<NewBoundaryRegionY>(mesh, false, mesh.iterateBndryUpperY()));
+      for (auto& bndry : mesh.getBoundaries()) {
+        if ((lower_y && bndry->location == BndryLoc::ydown)
+            or (upper_y && bndry->location == BndryLoc::yup)) {
+          boundary_regions.push_back(
+              std::dynamic_pointer_cast<bout::boundary::BoundaryRegionY>(bndry));
+        }
       }
     }
     // Cache boundary regions
     _contains.emplace_back(&mesh, false);
     _contains.emplace_back(&mesh, false);
-    iter_points([&](const auto& point) {
+    iter([&](const auto& point) {
       if (point.dir() == 1) {
         _contains[1][point.ind()] = true;
       } else if (point.dir() == -1) {
@@ -138,8 +127,17 @@ public:
   }
 
 private:
+  template <class Func>
+  void iter_regions(const Func& func) {
+    for (auto& region : boundary_regions) {
+      func(*region);
+    }
+    for (auto& region : boundary_regions_par) {
+      func(*region);
+    }
+  }
   std::vector<std::shared_ptr<bout::boundary::BoundaryRegionFCI>> boundary_regions_par;
-  std::vector<std::shared_ptr<NewBoundaryRegionY>> boundary_regions;
+  std::vector<std::shared_ptr<bout::boundary::BoundaryRegionY>> boundary_regions;
 
   std::vector<BoutMask> _contains;
 };
diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst
index 7ca782fc44..2aa464f7b7 100644
--- a/manual/sphinx/user_docs/boundary_options.rst
+++ b/manual/sphinx/user_docs/boundary_options.rst
@@ -503,8 +503,8 @@ geometries, as flux coordinate independent (FCI) method::
 
       void rhs() {
         BoutReal totalFlux = 0;
-        mesh->getCoordinates()->getYBoundary()->iter_points([&](auto& point) {
-          BoutReal flux = point.interpolate_sheath_o2(N) * point.interpolate_sheath_o2(V);
+        mesh->getCoordinates()->getYBoundary()->iter([&](auto& point) {
+          BoutReal flux = point.interpolate_boundary_o2(N) * point.interpolate_boundary_o2(V);
           totalFlux += flux;
         });
       }
@@ -531,18 +531,18 @@ Here is a short summary of some members of ``point``, where ``f`` is a :
 
    * - Function
      - Description
-   * - ``point.ythis(f)``
+   * - ``point.current(f)``
      - Returns the value at the last point in the domain
-   * - ``point.ynext(f)``
+   * - ``point.next(f)``
      - Returns the value at the first point in the boundary, i.e. one beyond the domain.
-   * - ``point.yprev(f)``
+   * - ``point.prev(f)``
      - Returns the value at the second to last point in the domain, if it is
        valid. NB: this point may not be valid.
-   * - ``point.interpolate_sheath_o2(f)``
+   * - ``point.interpolate_boundary_o2(f)``
      - Returns the value at the boundary, assuming the bounday value has been set
-   * - ``point.extrapolate_sheath_o1(f)``
+   * - ``point.extrapolate_boundary_o1(f)``
      - Returns the value at the boundary, extrapolating from the bulk, first order
-   * - ``point.extrapolate_sheath_o2(f)``
+   * - ``point.extrapolate_boundary_o2(f)``
      - Returns the value at the boundary, extrapolating from the bulk, second order
    * - ``point.extrapolate_next_o{1,2}(f)``
      - Extrapolate into the boundary from the bulk, first or second order
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index a513cc6021..8700096e24 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -619,7 +619,7 @@ void Field3D::applyParallelBoundary(const std::string& condition) {
   /// Loop over the mesh boundary regions
   for (const auto& reg : fieldmesh->getBoundariesPar()) {
     auto op = std::unique_ptr<BoundaryOpPar>{
-        dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg.get()))};
+        dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg))};
     op->apply(*this);
   }
 }
@@ -642,7 +642,7 @@ void Field3D::applyParallelBoundary(const std::string& region,
   for (const auto& reg : fieldmesh->getBoundariesPar()) {
     if (reg->label == region) {
       auto op = std::unique_ptr<BoundaryOpPar>{
-          dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg.get()))};
+          dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg))};
       op->apply(*this);
       break;
     }
@@ -669,7 +669,7 @@ void Field3D::applyParallelBoundary(const std::string& region,
       // BoundaryFactory can't create boundaries using Field3Ds, so get temporary
       // boundary of the right type
       auto tmp = std::unique_ptr<BoundaryOpPar>{
-          dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg.get()))};
+          dynamic_cast<BoundaryOpPar*>(bfact->create(condition, reg))};
       // then clone that with the actual argument
       auto op = std::unique_ptr<BoundaryOpPar>{tmp->clone(reg.get(), f)};
       op->apply(*this);
diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx
index 8f6a808482..f8935bb325 100644
--- a/src/field/field_data.cxx
+++ b/src/field/field_data.cxx
@@ -149,7 +149,7 @@ void FieldData::setBoundary(const std::string& name) {
   /// Get the mesh boundary regions
   /// Loop over the mesh parallel boundary regions
   for (const auto& reg : mesh->getBoundariesPar()) {
-    auto* op = dynamic_cast<BoundaryOpPar*>(bfact->createFromOptions(name, reg.get()));
+    auto* op = dynamic_cast<BoundaryOpPar*>(bfact->createFromOptions(name, reg));
     if (op != nullptr) {
       bndry_op_par.push_back(op);
     }
diff --git a/src/mesh/boundary_factory.cxx b/src/mesh/boundary_factory.cxx
index 988fcfce9b..f0c41cf8a1 100644
--- a/src/mesh/boundary_factory.cxx
+++ b/src/mesh/boundary_factory.cxx
@@ -11,6 +11,7 @@
 #include <array>
 #include <list>
 #include <map>
+#include <memory>
 #include <string>
 using std::list;
 using std::string;
@@ -85,7 +86,9 @@ void BoundaryFactory::cleanup() {
   instance = nullptr;
 }
 
-BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase* region) {
+BoundaryOpBase*
+BoundaryFactory::create(const string& name,
+                        const std::shared_ptr<BoundaryRegionBase>& region) {
 
   // Search for a string of the form: modifier(operation)
   auto pos = name.find('(');
@@ -107,8 +110,8 @@ BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase*
       // Clone the boundary operation, passing the region to operate over,
       // an empty args list and empty keyword map
       list<string> args;
-      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionFCI*>(region), args,
-                        {});
+      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionFCI*>(region.get()),
+                        args, {});
     } else {
       // Perpendicular boundary
       BoundaryOp* op = findBoundaryOp(trim(name));
@@ -199,16 +202,16 @@ BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase*
   if (pop != nullptr) {
     // An operation with arguments
     if (region->isParallel) {
-      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionFCI*>(region), arglist,
-                        keywords);
+      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionFCI*>(region.get()),
+                        arglist, keywords);
     }
     if (region->isX) {
-      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionX*>(region), arglist,
-                        keywords);
+      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionX*>(region.get()),
+                        arglist, keywords);
     }
     if (region->isY) {
-      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionY*>(region), arglist,
-                        keywords);
+      return pop->clone(dynamic_cast<bout::boundary::BoundaryRegionY*>(region.get()),
+                        arglist, keywords);
     }
   }
   if (!region->isParallel) {
@@ -227,12 +230,15 @@ BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase*
   return nullptr;
 }
 
-BoundaryOpBase* BoundaryFactory::create(const char* name, BoundaryRegionBase* region) {
+BoundaryOpBase*
+BoundaryFactory::create(const char* name,
+                        const std::shared_ptr<BoundaryRegionBase>& region) {
   return create(string(name), region);
 }
 
-BoundaryOpBase* BoundaryFactory::createFromOptions(const string& varname,
-                                                   BoundaryRegionBase* region) {
+BoundaryOpBase*
+BoundaryFactory::createFromOptions(const string& varname,
+                                   const std::shared_ptr<BoundaryRegionBase>& region) {
   if (region == nullptr) {
     return nullptr;
   }
@@ -357,8 +363,9 @@ BoundaryOpBase* BoundaryFactory::createFromOptions(const string& varname,
   // values. If a user want to override, specify "none" or "null"
 }
 
-BoundaryOpBase* BoundaryFactory::createFromOptions(const char* varname,
-                                                   BoundaryRegionBase* region) {
+BoundaryOpBase*
+BoundaryFactory::createFromOptions(const char* varname,
+                                   const std::shared_ptr<BoundaryRegionBase>& region) {
   return createFromOptions(string(varname), region);
 }
 
diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx
index 453e9c3d75..0e3c8bdcf3 100644
--- a/src/mesh/impls/bout/boutmesh.cxx
+++ b/src/mesh/impls/bout/boutmesh.cxx
@@ -94,11 +94,6 @@ BoutMesh::~BoutMesh() {
   // Delete the communication handles
   clear_handles();
 
-  // Delete the boundary regions
-  for (const auto& bndry : boundary) {
-    delete bndry;
-  }
-
   if (comm_x != MPI_COMM_NULL) {
     MPI_Comm_free(&comm_x);
   }
@@ -1177,7 +1172,9 @@ std::set<std::string> BoutMesh::getPossibleBoundaries() const {
         auto boundaries = mesh_copy.getBoundaries();
         std::transform(boundaries.begin(), boundaries.end(),
                        std::inserter(all_boundaries, all_boundaries.begin()),
-                       [](BoundaryRegionBase* boundary) { return boundary->label; });
+                       [](const std::shared_ptr<BoundaryRegionBase>& boundary) {
+                         return boundary->label;
+                       });
       };
 
   // This is sufficient to get the SOL boundary, if it exists
@@ -3249,7 +3246,9 @@ RangeIterator BoutMesh::iterateBndryUpperY() const {
   return RangeIterator(xs, xe);
 }
 
-std::vector<BoundaryRegionBase*> BoutMesh::getBoundaries() const { return boundary; }
+std::vector<std::shared_ptr<BoundaryRegionBase>> BoutMesh::getBoundaries() const {
+  return boundary;
+}
 
 using bout::boundary::BoundaryRegionFCI;
 std::vector<std::shared_ptr<BoundaryRegionFCI>>
diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx
index 79121ede58..a6d97a8944 100644
--- a/src/mesh/impls/bout/boutmesh.hxx
+++ b/src/mesh/impls/bout/boutmesh.hxx
@@ -169,7 +169,7 @@ public:
   bool hasBndryUpperY() const override { return has_boundary_upper_y; }
 
   // Boundary regions
-  std::vector<BoundaryRegionBase*> getBoundaries() const override;
+  std::vector<std::shared_ptr<BoundaryRegionBase>> getBoundaries() const override;
   std::vector<std::shared_ptr<bout::boundary::BoundaryRegionFCI>>
   getBoundariesPar(BoundaryParType type) const override;
   void addBoundaryPar(std::shared_ptr<bout::boundary::BoundaryRegionFCI> bndry,
@@ -416,7 +416,7 @@ protected:
   void addBoundaryRegions();
 
 private:
-  std::vector<BoundaryRegionBase*> boundary; // Vector of boundary regions
+  std::vector<std::shared_ptr<BoundaryRegionBase>> boundary; // Vector of boundary regions
   std::array<std::vector<std::shared_ptr<bout::boundary::BoundaryRegionFCI>>,
              static_cast<int>(BoundaryParType::SIZE)>
       par_boundary; // Vector of parallel boundary regions
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index e8c4da31ca..951230abc2 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -83,6 +83,7 @@ set(serial_tests_source
     ./include/bout/test_stencil.cxx
     ./include/bout/test_template_combinations.cxx
     ./include/bout/test_traits.cxx
+    ./include/bout/test_yboundary.cxx
     ./include/test_cyclic_reduction.cxx
     ./include/test_derivs.cxx
     ./include/test_mask.cxx
diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx
index 810e04cbb2..f957652a62 100644
--- a/tests/unit/fake_mesh.hxx
+++ b/tests/unit/fake_mesh.hxx
@@ -175,8 +175,12 @@ public:
   RangeIterator iterateBndryUpperInnerY() const override { return RangeIterator(); }
   bool hasBndryLowerY() const override { return false; }
   bool hasBndryUpperY() const override { return false; }
-  void addBoundary(BoundaryRegionBase* region) override { boundaries.push_back(region); }
-  std::vector<BoundaryRegionBase*> getBoundaries() const override { return boundaries; }
+  void addBoundary(std::shared_ptr<BoundaryRegionBase> region) override {
+    boundaries.push_back(region);
+  }
+  std::vector<std::shared_ptr<BoundaryRegionBase>> getBoundaries() const override {
+    return boundaries;
+  }
   std::vector<std::shared_ptr<bout::boundary::BoundaryRegionFCI>>
   getBoundariesPar(BoundaryParType UNUSED(type)) const override {
     return std::vector<std::shared_ptr<bout::boundary::BoundaryRegionFCI>>();
@@ -268,7 +272,7 @@ public:
   using Mesh::msg_len;
 
 private:
-  std::vector<BoundaryRegionBase*> boundaries;
+  std::vector<std::shared_ptr<BoundaryRegionBase>> boundaries;
 };
 
 /// FakeGridDataSource provides a non-null GridDataSource* source to use with FakeMesh, to
diff --git a/tests/unit/field/test_vector2d.cxx b/tests/unit/field/test_vector2d.cxx
index ded551084e..d99d2ed179 100644
--- a/tests/unit/field/test_vector2d.cxx
+++ b/tests/unit/field/test_vector2d.cxx
@@ -23,15 +23,8 @@ class Vector2DTest : public ::testing::Test {
 protected:
   Vector2DTest() {
     // Delete any existing mesh
-    if (mesh != nullptr) {
-      // Delete boundary regions
-      for (auto& r : mesh->getBoundaries()) {
-        delete r;
-      }
-
-      delete mesh;
-      mesh = nullptr;
-    }
+    delete mesh;
+    mesh = nullptr;
     bout::globals::mpi = new MpiWrapper();
     mesh = new FakeMesh(nx, ny, nz);
     static_cast<FakeMesh*>(mesh)->setCoordinates(nullptr);
@@ -60,14 +53,8 @@ class Vector2DTest : public ::testing::Test {
   }
 
   virtual ~Vector2DTest() {
-    if (mesh != nullptr) {
-      // Delete boundary regions
-      for (auto& r : mesh->getBoundaries()) {
-        delete r;
-      }
-      delete mesh;
-      mesh = nullptr;
-    }
+    delete mesh;
+    mesh = nullptr;
     delete mesh_staggered;
     mesh_staggered = nullptr;
     delete bout::globals::mpi;
diff --git a/tests/unit/field/test_vector3d.cxx b/tests/unit/field/test_vector3d.cxx
index d87d4c4bed..5b197d5cdf 100644
--- a/tests/unit/field/test_vector3d.cxx
+++ b/tests/unit/field/test_vector3d.cxx
@@ -21,11 +21,6 @@ class Vector3DTest : public ::testing::Test {
     WithQuietOutput quiet{output_info};
     // Delete any existing mesh
     if (mesh != nullptr) {
-      // Delete boundary regions
-      for (auto& r : mesh->getBoundaries()) {
-        delete r;
-      }
-
       delete mesh;
       mesh = nullptr;
     }
@@ -57,12 +52,6 @@ class Vector3DTest : public ::testing::Test {
   }
 
   ~Vector3DTest() override {
-    if (mesh != nullptr) {
-      // Delete boundary regions
-      for (auto& r : mesh->getBoundaries()) {
-        delete r;
-      }
-    }
     delete mesh;
     mesh = nullptr;
     delete mesh_staggered;
diff --git a/tests/unit/include/bout/test_yboundary.cxx b/tests/unit/include/bout/test_yboundary.cxx
new file mode 100644
index 0000000000..ce531c27f8
--- /dev/null
+++ b/tests/unit/include/bout/test_yboundary.cxx
@@ -0,0 +1,216 @@
+#include "gtest/gtest.h"
+
+#include "bout/boutexception.hxx"
+#include "bout/output_bout_types.hxx"
+#include "bout/yboundary_regions.hxx"
+
+#include "fake_mesh_fixture.hxx"
+#include "test_extras.hxx"
+
+using YBTest = FakeMeshFixture_tmpl<4, 5, 7>;
+
+using bout::globals::mesh;
+
+TEST_F(YBTest, dirichlet_o2_rgn) {
+  dynamic_cast<FakeMesh*>(mesh)->createBoundaries();
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { point.dirichlet_o2(test, 2); });
+  EXPECT_TRUE(IsFieldEqual(test, 3.0, "RGN_YGUARDS"));
+}
+
+TEST_F(YBTest, neumann_o1_rgn) {
+  dynamic_cast<FakeMesh*>(mesh)->createBoundaries();
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { point.neumann_o1(test, 1); });
+  EXPECT_TRUE(IsFieldEqual(test, 2.0, "RGN_YGUARDS"));
+}
+
+TEST_F(YBTest, neumann_o2_rgn) {
+  dynamic_cast<FakeMesh*>(mesh)->createBoundaries();
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { point.neumann_o2(test, 1); });
+  EXPECT_TRUE(IsFieldEqual(test, 3.0, "RGN_YGUARDS"));
+}
+
+TEST_F(YBTest, neumann_o3_rgn) {
+  dynamic_cast<FakeMesh*>(mesh)->createBoundaries();
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { point.neumann_o3(test, 1); });
+  EXPECT_TRUE(IsFieldEqual(test, 0.0, "RGN_YGUARDS"));
+}
+
+TEST_F(YBTest, bndry_size) {
+  dynamic_cast<FakeMesh*>(mesh)->createBoundaries();
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  int sum = 0;
+  sheath.iter([&](auto& point) { sum++; });
+  EXPECT_EQ(sum, 28);
+}
+
+TEST_F(YBTest, interpolate_boundary_o2) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter(
+      [&](auto& point) { EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 1.0); });
+}
+
+TEST_F(YBTest, interpolate_boundary_o2_const) {
+  const Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter(
+      [&](auto& point) { EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 1.0); });
+}
+
+TEST_F(YBTest, extrapolate_boundary_o2) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter(
+      [&](auto& point) { EXPECT_DOUBLE_EQ(point.extrapolate_boundary_o2(test), 1.0); });
+}
+
+TEST_F(YBTest, dirichlet_o1) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    point.dirichlet_o1(test, 2.0);
+    EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 1.5);
+  });
+}
+
+TEST_F(YBTest, dirichlet_o2) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    point.dirichlet_o2(test, 2.0);
+    EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 2.0);
+  });
+}
+
+TEST_F(YBTest, dirichlet_o3) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    point.dirichlet_o3(test, 2.0);
+    EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 7. / 3.);
+  });
+}
+
+TEST_F(YBTest, extrapolate_boundary_free) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2) + 1; }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    EXPECT_DOUBLE_EQ(point.extrapolate_boundary_free(
+                         test, bout::boundary::BoundaryFreeExtrapolation::limited),
+                     3);
+    EXPECT_DOUBLE_EQ(point.extrapolate_boundary_free(
+                         test, bout::boundary::BoundaryFreeExtrapolation::linear),
+                     3.5);
+    EXPECT_DOUBLE_EQ(point.extrapolate_boundary_free(
+                         test, bout::boundary::BoundaryFreeExtrapolation::exponential),
+                     5);
+  });
+}
+
+TEST_F(YBTest, set_free) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2) + 1; }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    point.set_free(test, bout::boundary::BoundaryFreeExtrapolation::limited);
+    EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 3);
+  });
+}
+
+TEST_F(YBTest, interpolate_boundary_o2_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter(
+      [&](auto& point) { EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 2.5); });
+}
+
+TEST_F(YBTest, limit_at_least) {
+  Field3D test = 1.0;
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    point.limit_at_least(test, 2.0);
+    EXPECT_DOUBLE_EQ(point.interpolate_boundary_o2(test), 1.5);
+  });
+}
+
+TEST_F(YBTest, extrapolate_grad_o2_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.extrapolate_grad_o2(test), 1); });
+}
+
+TEST_F(YBTest, extrapolate_next_o1_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.extrapolate_next_o1(test), 1); });
+}
+
+TEST_F(YBTest, extrapolate_next_o2_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.extrapolate_next_o2(test), 2); });
+}
+
+TEST_F(YBTest, next_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.next(test), 4); });
+}
+
+TEST_F(YBTest, current_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.current(test), 1); });
+}
+
+TEST_F(YBTest, prev_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.prev(test), 0); });
+}
+
+TEST_F(YBTest, next_const) {
+  const Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.next(test), 4); });
+}
+
+TEST_F(YBTest, current_const) {
+  const Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.current(test), 1); });
+}
+
+TEST_F(YBTest, prev_const) {
+  const Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) { EXPECT_DOUBLE_EQ(point.prev(test), 0); });
+}
+
+TEST_F(YBTest, getAt_square) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return SQ(i.y() - 2); }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  sheath.iter([&](auto& point) {
+    EXPECT_DOUBLE_EQ(point.getAt(test, 0), 4);
+    EXPECT_DOUBLE_EQ(point.getAt(test, 1), 1);
+    EXPECT_DOUBLE_EQ(point.getAt(test, 2), 0);
+  });
+}
+
+TEST_F(YBTest, getAt_func) {
+  Field3D test = makeField<Field3D>([&](auto& i) { return i.y() - 2; }, mesh);
+  YBoundary sheath(YBndryType::all, nullptr, *mesh);
+  auto square = [&](int yo, Ind3D ind) -> double { return SQ(test[ind]); };
+  sheath.iter([&](auto& point) {
+    EXPECT_DOUBLE_EQ(point.getAt(square, 0), 4);
+    EXPECT_DOUBLE_EQ(point.getAt(square, 1), 1);
+    EXPECT_DOUBLE_EQ(point.getAt(square, 2), 0);
+  });
+}
diff --git a/tests/unit/mesh/test_boundary_factory.cxx b/tests/unit/mesh/test_boundary_factory.cxx
index fd0100b755..fc9a6e95b0 100644
--- a/tests/unit/mesh/test_boundary_factory.cxx
+++ b/tests/unit/mesh/test_boundary_factory.cxx
@@ -39,21 +39,20 @@ class BoundaryFactoryTest : public ::testing::Test {
 
     fac->add(new TestBoundary(), "testboundary");
 
-    region = new BoundaryRegionXIn{"test_region", 0, 1, mesh};
+    region = std::make_shared<BoundaryRegionXIn>("test_region", 0, 1, mesh);
   }
 
   virtual ~BoundaryFactoryTest() {
     delete mesh;
     mesh = nullptr;
 
-    delete region;
     BoundaryFactory::cleanup();
 
     delete boundary;
   }
 
   BoundaryFactory* fac{BoundaryFactory::getInstance()};
-  BoundaryRegionXIn* region{nullptr};
+  std::shared_ptr<BoundaryRegionXIn> region;
   BoundaryOpBase* boundary{nullptr};
 };