diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp
index 14f966485ca5..6fdc42f572a1 100644
--- a/example/boundary_exchange/boundary_exchange.cpp
+++ b/example/boundary_exchange/boundary_exchange.cpp
@@ -12,9 +12,11 @@
 //========================================================================================
 
 // Standard Includes
+#include <cstdio>
 #include <iostream>
 #include <limits>
 #include <memory>
+#include <set>
 #include <string>
 #include <utility>
 #include <vector>
@@ -101,6 +103,185 @@ TaskStatus SetBlockValues(MeshData<Real> *md) {
   return TaskStatus::complete;
 }
 
+struct ParameterizedLine {
+  enum class ltype { straight, arc };
+
+  const Real x1, y1;
+  const Real x2, y2;
+  Real r, delta, phi;
+  ltype type;
+
+  using node = parthenon::forest::Node;
+  ParameterizedLine(std::shared_ptr<node> start, std::shared_ptr<node> end)
+      : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]},
+        type{ltype::straight} {
+    Real d1 = std::sqrt(x1 * x1 + y1 * y1);
+    Real d2 = std::sqrt(x2 * x2 + y2 * y2);
+    if (std::abs(d1 - d2) < 1.e-8 && d1 > 1.1) {
+      type = ltype::arc;
+      r = d1;
+      delta = M_PI / 4.0;
+      phi = 0.0;
+      if (y1 > 1.e-5) phi = delta;
+    }
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Real GetX(Real u) const {
+    if (type == ltype::straight)
+      return x1 * (1.0 - u) + x2 * u;
+    else
+      return r * cos(delta * u + phi);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Real GetY(Real u) const {
+    if (type == ltype::straight)
+      return y1 * (1.0 - u) + y2 * u;
+    else
+      return r * sin(delta * u + phi);
+  }
+};
+
+TaskStatus SetCoordinates(MeshData<Real> *md) {
+  using TE = parthenon::TopologicalElement;
+  auto pmesh = md->GetMeshPointer();
+  auto desc = parthenon::MakePackDescriptor<position>(md);
+
+  for (auto &ptree : pmesh->forest.GetTrees()) {
+    auto tree_id = ptree->GetId();
+    auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId()));
+    int i{0};
+    Real posx[4], posy[4];
+    for (auto &pnode : ptree->forest_nodes) {
+      posx[i] = pnode->x[0];
+      posy[i] = pnode->x[1];
+    }
+
+    auto &pnodes = ptree->forest_nodes;
+    ParameterizedLine c1(pnodes[0], pnodes[1]);
+    ParameterizedLine c3(pnodes[2], pnodes[3]);
+    ParameterizedLine c2(pnodes[0], pnodes[2]);
+    ParameterizedLine c4(pnodes[1], pnodes[3]);
+
+    IndexRange ib = md->GetBoundsI(IndexDomain::interior, TE::NN);
+    IndexRange jb = md->GetBoundsJ(IndexDomain::interior, TE::NN);
+    IndexRange kb = md->GetBoundsK(IndexDomain::interior, TE::NN);
+
+    parthenon::par_for(
+        parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0,
+        pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          auto &coords = pack.GetCoordinates(b);
+          const Real u = coords.X<X1DIR, TE::NN>(k, j, i);
+          const Real v = coords.X<X2DIR, TE::NN>(k, j, i);
+          const Real mu = 1.0 - u;
+          const Real mv = 1.0 - v;
+
+          Real x = mv * c1.GetX(u) + v * c3.GetX(u) + mu * c2.GetX(v) + u * c4.GetX(v) -
+                   (mu * mv * c1.GetX(0.0) + u * mv * c1.GetX(1.0) +
+                    mu * v * c2.GetX(1.0) + u * v * c3.GetX(1.0));
+          Real y = mv * c1.GetY(u) + v * c3.GetY(u) + mu * c2.GetY(v) + u * c4.GetY(v) -
+                   (mu * mv * c1.GetY(0.0) + u * mv * c1.GetY(1.0) +
+                    mu * v * c2.GetY(1.0) + u * v * c3.GetY(1.0));
+          pack(b, TE::NN, position(0), k, j, i) = x; // x-position
+          pack(b, TE::NN, position(1), k, j, i) = y; // y-position
+        });
+  }
+  return TaskStatus::complete;
+}
+
+TaskStatus FixTrivalentNodes2D(MeshData<Real> *md, bool coarse) {
+  using TE = parthenon::TopologicalElement;
+  auto pmesh = md->GetMeshPointer();
+  std::set<parthenon::PDOpt> opts{};
+  if (coarse) opts.insert(parthenon::PDOpt::Coarse);
+
+  auto desc = parthenon::MakePackDescriptor<position>(md, {}, opts);
+
+  const parthenon::CellLevel clevel =
+      coarse ? parthenon::CellLevel::coarse : parthenon::CellLevel::same;
+  IndexRange ib_in = md->GetBoundsI(clevel, IndexDomain::interior, TE::NN);
+  IndexRange jb_in = md->GetBoundsJ(clevel, IndexDomain::interior, TE::NN);
+  IndexRange kb = md->GetBoundsK(clevel, IndexDomain::interior, TE::NN);
+
+  for (auto &ptree : pmesh->forest.GetTrees()) {
+    auto tree_id = ptree->GetId();
+    int pos{0};
+    for (auto &pnode : ptree->forest_nodes) {
+      bool trivalent =
+          (pnode->associated_faces.size() == 3) && !pnode->on_physical_boundary;
+      if (trivalent) {
+        parthenon::CellCentOffsets offset(2 * (pos % 2) - 1, 2 * (pos / 2) - 1, 0);
+        auto pack =
+            desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId(), offset));
+        if (pack.GetNBlocks() == 0) return TaskStatus::complete;
+        // Copy shared boundary data
+
+        // This pack has to have only one block
+        auto gid = pack.GetGIDHost(0);
+        parthenon::MeshBlock *pmb;
+        for (auto &pmbd : md->GetAllBlockData())
+          if (pmbd->GetParentPointer()->gid == gid) pmb = pmbd->GetParentPointer();
+
+        // Now check which row of corner shared elements has been set by an owning block
+        parthenon::CellCentOffsets offsetX1(2 * (pos % 2) - 1, 0, 0);
+        parthenon::CellCentOffsets offsetX2(0, 2 * (pos / 2) - 1, 0);
+        int dir = X2DIR;
+        for (auto &neighbor : pmb->neighbors) {
+          if (neighbor.offsets == offsetX1 && neighbor.origin_ownership(offsetX2))
+            dir = X1DIR;
+          if (neighbor.offsets == offsetX2 && neighbor.origin_ownership(offsetX1))
+            dir = X2DIR;
+        }
+
+        // Select the reference location on the node shared by all three blocks
+        const int icorner = (offset(X1DIR) == parthenon::Offset::Low) ? ib_in.s : ib_in.e;
+        const int jcorner = (offset(X2DIR) == parthenon::Offset::Low) ? jb_in.s : jb_in.e;
+        // Select the index space of elements that need to get overwritten
+        IndexRange ib = offset(X1DIR) == parthenon::Offset::Low
+                            ? IndexRange{0, ib_in.s - (dir == X2DIR)}
+                            : IndexRange{ib_in.e + (dir == X2DIR),
+                                         ib_in.e + parthenon::Globals::nghost};
+        IndexRange jb = offset(X2DIR) == parthenon::Offset::Low
+                            ? IndexRange{0, jb_in.s - (dir == X1DIR)}
+                            : IndexRange{jb_in.e + (dir == X1DIR),
+                                         ib_in.e + parthenon::Globals::nghost};
+        parthenon::par_for(
+            parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0,
+            pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+            KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+              const int enclosed_area = std::abs((i - icorner) * (j - jcorner));
+              int ioff{0};
+              int joff{0};
+              // This is an element shared by both of the neighbors, but only one
+              // communicated it, so copy it to the other side of the corner
+              if (enclosed_area == 0) {
+                joff = std::abs(icorner - i) * offset(X2DIR);
+                ioff = std::abs(jcorner - j) * offset(X1DIR);
+              } else {
+                // These are non-existent elements that we fill with data copied from
+                // their nearest neighbors
+                int off = 1;
+                while (off * off < enclosed_area)
+                  off++;
+                ioff = offset(X1DIR) * off * (dir == X1DIR);
+                joff = offset(X2DIR) * off * (dir == X2DIR);
+              }
+              pack(b, TE::NN, position(0), k, j, i) =
+                  pack(b, TE::NN, position(0), k, jcorner + joff,
+                       icorner + ioff); // x-position
+              pack(b, TE::NN, position(1), k, j, i) =
+                  pack(b, TE::NN, position(1), k, jcorner + joff,
+                       icorner + ioff); // y-position
+            });
+      }
+      pos++;
+    }
+  }
+  return TaskStatus::complete;
+}
+
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto package = std::make_shared<StateDescriptor>("boundary_exchange");
   Params &params = package->AllParams();
@@ -109,8 +290,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
              std::vector<int>{8});
   m.RegisterRefinementOps<parthenon::refinement_ops::ProlongatePiecewiseConstant,
                           parthenon::refinement_ops::RestrictAverage>();
-  package->AddField(neighbor_info::name(), m);
+  package->AddField<neighbor_info>(m);
 
+  Metadata m_node({Metadata::Node, Metadata::Independent, Metadata::FillGhost},
+                  std::vector<int>{2});
+  package->AddField<position>(m_node);
   return package;
 }
 
diff --git a/example/boundary_exchange/boundary_exchange.hpp b/example/boundary_exchange/boundary_exchange.hpp
index 7b0cbb87bd88..bd8a60cbfcdc 100644
--- a/example/boundary_exchange/boundary_exchange.hpp
+++ b/example/boundary_exchange/boundary_exchange.hpp
@@ -34,7 +34,16 @@ struct neighbor_info : public parthenon::variable_names::base_t<false, 4> {
   static std::string name() { return "neighbor_info"; }
 };
 
+struct position : public parthenon::variable_names::base_t<false, 4> {
+  template <class... Ts>
+  KOKKOS_INLINE_FUNCTION position(Ts &&...args)
+      : parthenon::variable_names::base_t<false, 4>(std::forward<Ts>(args)...) {}
+  static std::string name() { return "position"; }
+};
+
 TaskStatus SetBlockValues(MeshData<Real> *rc);
+TaskStatus SetCoordinates(MeshData<Real> *rc);
+TaskStatus FixTrivalentNodes2D(MeshData<Real> *rc, bool coarse);
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
 
 } // namespace boundary_exchange
diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp
index 640705ca77ef..0b0fc408a7d4 100644
--- a/example/boundary_exchange/boundary_exchange_driver.cpp
+++ b/example/boundary_exchange/boundary_exchange_driver.cpp
@@ -34,24 +34,7 @@ using boundary_exchange::BoundaryExchangeDriver;
 
 Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin);
 
-int main(int argc, char *argv[]) {
-  ParthenonManager pman;
-
-  pman.app_input->ProcessPackages = ProcessPackages;
-
-  // This is called on each mesh block whenever the mesh changes.
-  // pman.app_input->InitMeshBlockUserData = &calculate_pi::SetInOrOutBlock;
-
-  auto manager_status = pman.ParthenonInitEnv(argc, argv);
-  if (manager_status == ParthenonStatus::complete) {
-    pman.ParthenonFinalize();
-    return 0;
-  }
-  if (manager_status == ParthenonStatus::error) {
-    pman.ParthenonFinalize();
-    return 1;
-  }
-
+parthenon::forest::ForestDefinition MakeFourTreeForestWithOneTreeRotated() {
   // Create the nodes for the forest, the x-y positions are only used for
   // visualizing the forest configuration and *do not* determine the global
   // coordinates of the trees.
@@ -95,7 +78,100 @@ int main(int argc, char *argv[]) {
   forest_def.AddBC(edge_t({n[7], n[8]}));
 
   forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, 1, 0, 0, 0));
-  pman.ParthenonInitPackagesAndMesh(forest_def);
+
+  return forest_def;
+}
+
+parthenon::forest::ForestDefinition n_blocks(int nblocks) {
+  using namespace parthenon::forest;
+  std::unordered_map<uint64_t, std::shared_ptr<Node>> nodes;
+  ForestDefinition forest_def;
+  parthenon::Real xoffset = 0.1;
+  parthenon::Real yoffset = 0.1;
+
+  for (int point = 0; point < 2 * nblocks; ++point) {
+    nodes[point] = Node::create(point, {-std::sin(point * M_PI / nblocks) + xoffset,
+                                        -std::cos(point * M_PI / nblocks) + yoffset});
+  }
+
+  using edge_t = parthenon::forest::Edge;
+  for (int point = 0; point < 2 * nblocks; ++point) {
+    forest_def.AddBC(
+        edge_t({nodes[point % (2 * nblocks)], nodes[(point + 1) % (2 * nblocks)]}));
+  }
+
+  nodes[2 * nblocks] = Node::create(2 * nblocks, {xoffset, yoffset});
+  auto &n = nodes;
+  for (int t = 0; t < nblocks; ++t)
+    forest_def.AddFace(
+        t, {n[2 * t + 1], n[2 * t], n[(2 * t + 2) % (2 * nblocks)], n[2 * nblocks]});
+
+  int level = 1;
+  forest_def.AddInitialRefinement(
+      parthenon::LogicalLocation(0, level, (1 << level) - 1, (1 << level) - 1, 0));
+  forest_def.AddInitialRefinement(
+      parthenon::LogicalLocation(1, level, (1 << level) - 1, (1 << level) - 1, 0));
+  forest_def.AddInitialRefinement(
+      parthenon::LogicalLocation(2, level, (1 << level) - 1, (1 << level) - 1, 0));
+
+  return forest_def;
+}
+
+parthenon::forest::ForestDefinition cubed_circle(std::vector<double> radii) {
+  using namespace parthenon::forest;
+  std::vector<std::shared_ptr<Node>> nodes_x, nodes_y, nodes_c;
+  std::size_t node_idx{0};
+  std::shared_ptr<Node> origin = Node::create(node_idx++, {0.0, 0.0});
+  double fac = 1.0 / sqrt(2.0);
+  for (auto rad : radii) {
+    nodes_x.push_back(Node::create(node_idx++, {rad, 0.0}));
+    nodes_c.push_back(Node::create(node_idx++, {rad * fac, rad * fac}));
+    nodes_y.push_back(Node::create(node_idx++, {0.0, rad}));
+  }
+
+  ForestDefinition forest_def;
+  std::size_t face_idx{0};
+  forest_def.AddFace(face_idx++, {origin, nodes_x[0], nodes_y[0], nodes_c[0]});
+  forest_def.AddBC(Edge({origin, nodes_x[0]}));
+  forest_def.AddBC(Edge({origin, nodes_y[0]}));
+  for (int i = 0; i < nodes_x.size() - 1; ++i) {
+    forest_def.AddFace(face_idx++,
+                       {nodes_x[i], nodes_x[i + 1], nodes_c[i], nodes_c[i + 1]});
+    forest_def.AddFace(face_idx++,
+                       {nodes_c[i], nodes_c[i + 1], nodes_y[i], nodes_y[i + 1]});
+    forest_def.AddBC(Edge({nodes_x[i], nodes_x[i + 1]}));
+    forest_def.AddBC(Edge({nodes_y[i], nodes_y[i + 1]}));
+  }
+  forest_def.AddBC(Edge({nodes_x.back(), nodes_c.back()}));
+  forest_def.AddBC(Edge({nodes_c.back(), nodes_y.back()}));
+
+  int level = 2;
+  forest_def.AddInitialRefinement(
+      parthenon::LogicalLocation(1, level, 0, (1 << level) - 1, 0));
+
+  return forest_def;
+}
+
+int main(int argc, char *argv[]) {
+  ParthenonManager pman;
+
+  pman.app_input->ProcessPackages = ProcessPackages;
+
+  // This is called on each mesh block whenever the mesh changes.
+  // pman.app_input->InitMeshBlockUserData = &calculate_pi::SetInOrOutBlock;
+
+  auto manager_status = pman.ParthenonInitEnv(argc, argv);
+  if (manager_status == ParthenonStatus::complete) {
+    pman.ParthenonFinalize();
+    return 0;
+  }
+  if (manager_status == ParthenonStatus::error) {
+    pman.ParthenonFinalize();
+    return 1;
+  }
+
+  pman.ParthenonInitPackagesAndMesh(cubed_circle({1.0, 2.0, 3.0}));
+  // pman.ParthenonInitPackagesAndMesh(n_blocks(3));
 
   // This needs to be scoped so that the driver object is destructed before Finalize
   {
@@ -150,7 +226,14 @@ TaskCollection BoundaryExchangeDriver::MakeTaskCollection(T &blocks) {
       auto &md = pmesh->mesh_data.GetOrAdd("base", i);
       TaskID none(0);
       auto fill = tl.AddTask(none, SetBlockValues, md.get());
-      auto bound = AddBoundaryExchangeTasks(fill, tl, md, true);
+      auto set_coords = tl.AddTask(none, SetCoordinates, md.get());
+      auto bound = AddBoundaryExchangeTasks(
+          fill | set_coords, tl, md, true,
+          [](TaskID depend, TaskList *ptl, std::shared_ptr<MeshData<Real>> md,
+             bool coarse) {
+            return ptl->AddTask(depend, FixTrivalentNodes2D, md.get(), coarse);
+          });
+      // auto fix = tl.AddTask(bound, FixTrivalentNodes2D, md.get(), false);
     }
   }
 
diff --git a/src/bvals/comms/bnd_info.cpp b/src/bvals/comms/bnd_info.cpp
index 2b4948278c38..370dcac9962b 100644
--- a/src/bvals/comms/bnd_info.cpp
+++ b/src/bvals/comms/bnd_info.cpp
@@ -244,7 +244,12 @@ CalcIndices(const NeighborBlock &nb, MeshBlock *pmb,
       if (sox2 == 0) sox2 = loc.l(1) % 2 == 1 ? 1 : -1;
       if (sox3 == 0) sox3 = loc.l(2) % 2 == 1 ? 1 : -1;
     }
-    owns = GetIndexRangeMaskFromOwnership(el, nb.ownership, sox1, sox2, sox3);
+    // We are working in the frame of the origin block (i.e. the receiving block),
+    // so we use the ownership array in that frame to get the masking initially.
+    // This mask must then be transformed to the frame of the neighbor block (i.e.
+    // sending block) since the index range we are masking is defined in that frame.
+    owns = lcoord_trans.Transform(
+        GetIndexRangeMaskFromOwnership(el, nb.origin_ownership, sox1, sox2, sox3));
   }
   return SpatiallyMaskedIndexer6D(owns, {0, tensor_shape[0] - 1},
                                   {0, tensor_shape[1] - 1}, {0, tensor_shape[2] - 1},
diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp
index 78121cd3fac9..222c1f5c58fc 100644
--- a/src/bvals/comms/boundary_communication.cpp
+++ b/src/bvals/comms/boundary_communication.cpp
@@ -302,10 +302,11 @@ TaskStatus SetBounds(std::shared_ptr<MeshData<Real>> &md) {
                   const int ii = i;
                   Kokkos::parallel_for(
                       Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) {
-                        const auto [il, jl, kl] =
-                            lcoord_trans.InverseTransform({ii + m, jj, kk});
-                        if (idxer.IsActive(kl, jl, il))
+                        if (idxer.IsActive(kk, jj, ii + m)) {
+                          const auto [il, jl, kl] = lcoord_trans.InverseTransform(
+                              std::array<int, 3>{ii + m, jj, kk});
                           var(iel, tt, uu, vv, kl, jl, il) = fac * buf[m];
+                        }
                       });
                 });
           } else if (bnd_info(b).allocated && bound_type != BoundaryType::flxcor_recv) {
@@ -322,10 +323,11 @@ TaskStatus SetBounds(std::shared_ptr<MeshData<Real>> &md) {
                   const int ii = i;
                   Kokkos::parallel_for(
                       Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) {
-                        const auto [il, jl, kl] =
-                            lcoord_trans.InverseTransform({ii + m, jj, kk});
-                        if (idxer.IsActive(kl, jl, il))
+                        if (idxer.IsActive(kk, jj, ii + m)) {
+                          const auto [il, jl, kl] = lcoord_trans.InverseTransform(
+                              std::array<int, 3>{ii + m, jj, kk});
                           var(iel, tt, uu, vv, kl, jl, il) = default_val;
+                        }
                       });
                 });
           }
@@ -401,10 +403,77 @@ ProlongateBounds<BoundaryType::nonlocal>(std::shared_ptr<MeshData<Real>> &);
 template TaskStatus
 ProlongateBounds<BoundaryType::gmg_prolongate_recv>(std::shared_ptr<MeshData<Real>> &);
 
+template <BoundaryType bound_type>
+TaskStatus ProlongateInternalBounds(std::shared_ptr<MeshData<Real>> &md) {
+  PARTHENON_INSTRUMENT
+
+  Mesh *pmesh = md->GetMeshPointer();
+  auto &cache = md->GetBvarsCache().GetSubCache(bound_type, false);
+
+  auto [rebuild, nbound] = CheckReceiveBufferCacheForRebuild<bound_type, false>(md);
+
+  if (rebuild) {
+    if constexpr (bound_type == BoundaryType::gmg_prolongate_recv) {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetInteriorProlongate);
+    } else if constexpr (bound_type == BoundaryType::gmg_restrict_recv) {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetNull);
+    } else {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetSet);
+    }
+  }
+
+  if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) {
+    auto pmb = md->GetBlockData(0)->GetBlockPointer();
+    StateDescriptor *resolved_packages = pmb->resolved_packages.get();
+
+    // Prolongate from coarse buffer
+    refinement::ProlongateInternal(resolved_packages, cache.prores_cache, pmb->cellbounds,
+                                   pmb->c_cellbounds);
+  }
+  return TaskStatus::complete;
+}
+
+template <BoundaryType bound_type>
+TaskStatus ProlongateSharedBounds(std::shared_ptr<MeshData<Real>> &md) {
+  PARTHENON_INSTRUMENT
+
+  Mesh *pmesh = md->GetMeshPointer();
+  auto &cache = md->GetBvarsCache().GetSubCache(bound_type, false);
+
+  auto [rebuild, nbound] = CheckReceiveBufferCacheForRebuild<bound_type, false>(md);
+
+  if (rebuild) {
+    if constexpr (bound_type == BoundaryType::gmg_prolongate_recv) {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetInteriorProlongate);
+    } else if constexpr (bound_type == BoundaryType::gmg_restrict_recv) {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetNull);
+    } else {
+      RebuildBufferCache<bound_type, false>(md, nbound, BndInfo::GetSetBndInfo,
+                                            ProResInfo::GetSet);
+    }
+  }
+
+  if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) {
+    auto pmb = md->GetBlockData(0)->GetBlockPointer();
+    StateDescriptor *resolved_packages = pmb->resolved_packages.get();
+
+    // Prolongate from coarse buffer
+    refinement::ProlongateShared(resolved_packages, cache.prores_cache, pmb->cellbounds,
+                                 pmb->c_cellbounds);
+  }
+  return TaskStatus::complete;
+}
+
 // Adds all relevant boundary communication to a single task list
 template <BoundaryType bounds>
 TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl,
-                                std::shared_ptr<MeshData<Real>> &md, bool multilevel) {
+                                std::shared_ptr<MeshData<Real>> &md, bool multilevel,
+                                bound_task_adder_t extra_bounds) {
   // TODO(LFR): Splitting up the boundary tasks while doing prolongation can cause some
   //            possible issues for sparse fields. In particular, the order in which
   //            fields are allocated and then set could potentially result in different
@@ -437,19 +506,22 @@ TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl,
   auto pro = set;
   if (md->GetMeshPointer()->multilevel) {
     auto cbound = tl.AddTask(set, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, true);
+    if (extra_bounds) cbound = extra_bounds(cbound, &tl, md, true);
     pro = tl.AddTask(cbound, TF(ProlongateBounds<bounds>), md);
+    auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false);
+    if (extra_bounds) fbound = extra_bounds(fbound, &tl, md, false);
+    return tl.AddTask(fbound, TF(ProlongateInternalBounds<bounds>), md);
+  } else {
+    auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false);
+    if (extra_bounds) fbound = extra_bounds(fbound, &tl, md, false);
+    return fbound;
   }
-  auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false);
-
-  return fbound;
 }
-template TaskID
-AddBoundaryExchangeTasks<BoundaryType::any>(TaskID, TaskList &,
-                                            std::shared_ptr<MeshData<Real>> &, bool);
+template TaskID AddBoundaryExchangeTasks<BoundaryType::any>(
+    TaskID, TaskList &, std::shared_ptr<MeshData<Real>> &, bool, bound_task_adder_t);
 
-template TaskID
-AddBoundaryExchangeTasks<BoundaryType::gmg_same>(TaskID, TaskList &,
-                                                 std::shared_ptr<MeshData<Real>> &, bool);
+template TaskID AddBoundaryExchangeTasks<BoundaryType::gmg_same>(
+    TaskID, TaskList &, std::shared_ptr<MeshData<Real>> &, bool, bound_task_adder_t);
 
 TaskID AddFluxCorrectionTasks(TaskID dependency, TaskList &tl,
                               std::shared_ptr<MeshData<Real>> &md, bool multilevel) {
diff --git a/src/bvals/comms/bvals_in_one.hpp b/src/bvals/comms/bvals_in_one.hpp
index 0180bd6a8305..a38952ccc58e 100644
--- a/src/bvals/comms/bvals_in_one.hpp
+++ b/src/bvals/comms/bvals_in_one.hpp
@@ -80,9 +80,12 @@ static TaskStatus SetFluxCorrections(std::shared_ptr<MeshData<Real>> &md) {
 }
 
 // Adds all relevant boundary communication to a single task list
+using bound_task_adder_t =
+    std::function<TaskID(TaskID, TaskList *, std::shared_ptr<MeshData<Real>>, bool)>;
 template <BoundaryType bounds = BoundaryType::any>
 TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl,
-                                std::shared_ptr<MeshData<Real>> &md, bool multilevel);
+                                std::shared_ptr<MeshData<Real>> &md, bool multilevel,
+                                bound_task_adder_t extra_bounds = bound_task_adder_t());
 
 // Adds all relevant flux correction tasks to a single task list
 TaskID AddFluxCorrectionTasks(TaskID dependency, TaskList &tl,
diff --git a/src/bvals/neighbor_block.hpp b/src/bvals/neighbor_block.hpp
index 5dfc235ff3a5..da96a910bce9 100644
--- a/src/bvals/neighbor_block.hpp
+++ b/src/bvals/neighbor_block.hpp
@@ -63,7 +63,10 @@ struct NeighborBlock {
   // Offset of the neighbor block relative to origin block
   CellCentOffsets offsets;
   // Ownership of neighbor block of different topological elements
+  // (in the frame of the neighbor block)
   block_ownership_t ownership;
+  // (in the frame of the origin block)
+  block_ownership_t origin_ownership;
   // Logical coordinate transformation from the main block to this neighbor
   forest::LogicalCoordinateTransformation lcoord_trans;
 
diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp
index db9de2e0025b..7f4d93931ae9 100644
--- a/src/mesh/forest/forest.hpp
+++ b/src/mesh/forest/forest.hpp
@@ -68,6 +68,8 @@ class ForestDefinition {
       PARTHENON_REQUIRE(periodic_connection,
                         "Must specify another edge for periodic boundary conditions.");
     bc_edges.emplace_back(ForestBC<Edge>{edge, bf, periodic_connection});
+    for (auto &pnode : edge.nodes)
+      pnode->on_physical_boundary = true;
   }
 
   void AddInitialRefinement(const LogicalLocation &loc) {
diff --git a/src/mesh/forest/forest_node.hpp b/src/mesh/forest/forest_node.hpp
index d95e2a59ee89..0d523f63ea98 100644
--- a/src/mesh/forest/forest_node.hpp
+++ b/src/mesh/forest/forest_node.hpp
@@ -33,7 +33,8 @@ namespace forest {
 class Face;
 class Node {
  public:
-  Node(int id_in, std::array<Real, NDIM> pos) : id(id_in), x(pos), associated_faces{} {}
+  Node(int id_in, std::array<Real, NDIM> pos)
+      : id(id_in), x(pos), associated_faces{}, on_physical_boundary{false} {}
 
   static std::shared_ptr<Node> create(int id, std::array<Real, NDIM> pos) {
     return std::make_shared<Node>(id, pos);
@@ -43,6 +44,7 @@ class Node {
   std::array<Real, NDIM> x;
   std::unordered_set<std::weak_ptr<Face>, WeakPtrHash<Face>, WeakPtrEqual<Face>>
       associated_faces;
+  bool on_physical_boundary;
 };
 } // namespace forest
 } // namespace parthenon
diff --git a/src/mesh/forest/logical_coordinate_transformation.cpp b/src/mesh/forest/logical_coordinate_transformation.cpp
index 34d4bec27a6f..55242844c867 100644
--- a/src/mesh/forest/logical_coordinate_transformation.cpp
+++ b/src/mesh/forest/logical_coordinate_transformation.cpp
@@ -100,6 +100,44 @@ CellCentOffsets LogicalCoordinateTransformation::Transform(CellCentOffsets in) c
   return out;
 }
 
+CellCentOffsets
+LogicalCoordinateTransformation::InverseTransform(CellCentOffsets in) const {
+  CellCentOffsets out;
+  for (int dir = 0; dir < 3; ++dir) {
+    const int indir = abs(dir_connection[dir]);
+    out.u[dir] = dir_flip[dir] ? -in.u[indir] : in.u[indir];
+  }
+  return out;
+}
+
+block_ownership_t
+LogicalCoordinateTransformation::Transform(const block_ownership_t &in) const {
+  block_ownership_t out;
+  for (int ox3 : {-1, 0, 1}) {
+    for (int ox2 : {-1, 0, 1}) {
+      for (int ox1 : {-1, 0, 1}) {
+        auto offset = CellCentOffsets(ox1, ox2, ox3);
+        out(Transform(offset)) = in(offset);
+      }
+    }
+  }
+  return out;
+}
+
+block_ownership_t
+LogicalCoordinateTransformation::InverseTransform(const block_ownership_t &in) const {
+  block_ownership_t out;
+  for (int ox3 : {-1, 0, 1}) {
+    for (int ox2 : {-1, 0, 1}) {
+      for (int ox1 : {-1, 0, 1}) {
+        auto offset = CellCentOffsets(ox1, ox2, ox3);
+        out(offset) = in(Transform(offset));
+      }
+    }
+  }
+  return out;
+}
+
 LogicalCoordinateTransformation
 ComposeTransformations(const LogicalCoordinateTransformation &first,
                        const LogicalCoordinateTransformation &second) {
diff --git a/src/mesh/forest/logical_coordinate_transformation.hpp b/src/mesh/forest/logical_coordinate_transformation.hpp
index 0e0d0bc293cf..57aec7db9691 100644
--- a/src/mesh/forest/logical_coordinate_transformation.hpp
+++ b/src/mesh/forest/logical_coordinate_transformation.hpp
@@ -51,6 +51,9 @@ struct LogicalCoordinateTransformation {
   LogicalLocation InverseTransform(const LogicalLocation &loc_in,
                                    std::int64_t origin) const;
   CellCentOffsets Transform(CellCentOffsets in) const;
+  CellCentOffsets InverseTransform(CellCentOffsets in) const;
+  block_ownership_t Transform(const block_ownership_t &in) const;
+  block_ownership_t InverseTransform(const block_ownership_t &in) const;
 
   KOKKOS_INLINE_FUNCTION
   std::tuple<TopologicalElement, Real> Transform(TopologicalElement el) const {
@@ -112,9 +115,12 @@ struct NeighborLocation {
   NeighborLocation(const LogicalLocation &g, const LogicalLocation &o,
                    const LogicalCoordinateTransformation &lcoord_trans)
       : global_loc(g), origin_loc(o), lcoord_trans(lcoord_trans) {}
-  LogicalLocation global_loc; // Global location of neighboring block
-  LogicalLocation
-      origin_loc; // Logical location of neighboring block in index space of origin block
+  // Global location of neighboring block
+  LogicalLocation global_loc;
+  // Logical location of neighboring block in index space of origin block
+  LogicalLocation origin_loc;
+  // Coordinate transform to coords of neighbor tree from origin tree
+  // (i.e. global_loc = lcoord_trans.Transform(origin_loc))
   LogicalCoordinateTransformation lcoord_trans;
 };
 
diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp
index 85afb715202d..8d1d3c5ee070 100644
--- a/src/mesh/forest/logical_location.hpp
+++ b/src/mesh/forest/logical_location.hpp
@@ -34,6 +34,7 @@
 
 #include "basic_types.hpp"
 #include "defs.hpp"
+#include "utils/cell_center_offsets.hpp"
 #include "utils/error_checking.hpp"
 #include "utils/morton_number.hpp"
 
@@ -123,6 +124,10 @@ class LogicalLocation { // aggregate and POD type
     return {(lx1() & 1LL) == 1LL, (lx2() & 1LL) == 1LL, (lx3() & 1LL) == 1LL};
   }
 
+  bool IsOnTreeBoundary(CellCentOffsets offset) const {
+    return IsOnTreeBoundary(offset(X1DIR), offset(X2DIR), offset(X3DIR));
+  }
+
   bool IsOnTreeBoundary(int ox1, int ox2, int ox3) const {
     const int nup = (1 << std::max(level(), 0)) - 1;
     const bool bound1 =
diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp
index e51ad8e5daec..f51747b057a2 100644
--- a/src/mesh/mesh-gmg.cpp
+++ b/src/mesh/mesh-gmg.cpp
@@ -76,13 +76,15 @@ void Mesh::SetMeshBlockNeighbors(
       all_neighbors.emplace_back(pmb->pmy_mesh, nloc.global_loc, nloc.origin_loc,
                                  ranklist[lgid], gid, offsets, bid, tid, f[0], f[1]);
 
-      // Set neighbor block ownership
+      // Set neighbor block ownership (note that this is ownership in coordinates of the
+      // neighbor block)
       auto &nb = all_neighbors.back();
       auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc, grid_id);
-
       nb.ownership =
           DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined);
       nb.ownership.initialized = true;
+      nb.origin_ownership = nloc.lcoord_trans.InverseTransform(nb.ownership);
+      nb.origin_ownership.initialized = true;
 
       // Set logical coordinate transformation from this block to the neighbor
       nb.lcoord_trans = nloc.lcoord_trans;
diff --git a/src/pack/block_selector.hpp b/src/pack/block_selector.hpp
index 1283ecffd573..71b25f1cf59a 100644
--- a/src/pack/block_selector.hpp
+++ b/src/pack/block_selector.hpp
@@ -32,6 +32,7 @@
 #include "interface/variable.hpp"
 #include "pack/pack_utils.hpp"
 #include "pack/sparse_pack_base.hpp"
+#include "utils/cell_center_offsets.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/type_list.hpp"
 #include "utils/utils.hpp"
@@ -74,6 +75,14 @@ inline block_selector_func_t OnPhysicalBoundary() {
   };
 }
 
+inline block_selector_func_t OnTree(std::size_t tree_id,
+                                    CellCentOffsets offset = CellCentOffsets(0, 0, 0)) {
+  return [tree_id, offset](MeshBlockData<Real> *pmbd) {
+    const auto &loc = pmbd->GetBlockPointer()->loc;
+    return (loc.tree() == tree_id) && loc.IsOnTreeBoundary(offset);
+  };
+}
+
 inline block_selector_func_t FineOnCompositeGrid(const MeshData<Real> *pmd) {
   if (pmd->grid.type == GridType::two_level_composite) {
     const int fine_level = pmd->grid.logical_level;
diff --git a/src/utils/cell_center_offsets.hpp b/src/utils/cell_center_offsets.hpp
index 1b27da0d9779..49cbbb651d6e 100644
--- a/src/utils/cell_center_offsets.hpp
+++ b/src/utils/cell_center_offsets.hpp
@@ -26,7 +26,6 @@
 #include "basic_types.hpp"
 #include "defs.hpp"
 #include "utils/bit_hacks.hpp"
-#include "utils/indexer.hpp"
 
 namespace parthenon {
 // CellCentOffsets defines the position of a topological element
@@ -39,7 +38,7 @@ namespace parthenon {
 
 // TODO(LFR): Consider switching this to C-style enum within a namespace to avoid
 // static_cast
-enum class Offset : int { Low = -1, Middle = 0, Up = 1 };
+enum Offset : int { Low = -1, Middle = 0, Up = 1 };
 inline int operator+(Offset a, int b) { return static_cast<int>(a) + b; }
 inline int operator+(int b, Offset a) { return static_cast<int>(a) + b; }
 inline Offset operator-(Offset in) { return static_cast<Offset>(-static_cast<int>(in)); }
@@ -58,7 +57,7 @@ struct CellCentOffsets {
 
   Offset &operator[](int idx) { return u[idx]; }
   const Offset &operator[](int idx) const { return u[idx]; }
-  int operator()(CoordinateDirection dir) const { return static_cast<int>(u[dir - 1]); }
+  Offset operator()(CoordinateDirection dir) const { return u[dir - 1]; }
 
   operator std::array<int, 3>() const {
     return {static_cast<int>(u[0]), static_cast<int>(u[1]), static_cast<int>(u[2])};
@@ -75,6 +74,8 @@ struct CellCentOffsets {
   // element in that direction from the cell center.
   std::vector<std::pair<CoordinateDirection, Offset>> GetNormals() const;
 
+  bool operator==(const CellCentOffsets &other) const { return u == other.u; }
+
   bool IsNode() const {
     return 3 == abs(static_cast<int>(u[0])) + abs(static_cast<int>(u[1])) +
                     abs(static_cast<int>(u[2]));
diff --git a/src/utils/indexer.hpp b/src/utils/indexer.hpp
index 723c090ac25f..e764e465628c 100644
--- a/src/utils/indexer.hpp
+++ b/src/utils/indexer.hpp
@@ -21,6 +21,7 @@
 
 #include <Kokkos_Core.hpp>
 
+#include "utils/cell_center_offsets.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/type_list.hpp"
 #include "utils/utils.hpp"
@@ -37,6 +38,14 @@ struct block_ownership_t {
   bool &operator()(int ox1, int ox2, int ox3) {
     return ownership[ox1 + 1][ox2 + 1][ox3 + 1];
   }
+  KOKKOS_FORCEINLINE_FUNCTION
+  bool &operator()(const CellCentOffsets &offset) {
+    return ownership[offset[0] + 1][offset[1] + 1][offset[2] + 1];
+  }
+  KOKKOS_FORCEINLINE_FUNCTION
+  const bool &operator()(const CellCentOffsets &offset) const {
+    return ownership[offset[0] + 1][offset[1] + 1][offset[2] + 1];
+  }
 
   KOKKOS_FORCEINLINE_FUNCTION
   block_ownership_t() : block_ownership_t(false) {}