diff --git a/src/atlas/functionspace/CellColumns.cc b/src/atlas/functionspace/CellColumns.cc
index 01294f6c5..7a65543c2 100644
--- a/src/atlas/functionspace/CellColumns.cc
+++ b/src/atlas/functionspace/CellColumns.cc
@@ -16,6 +16,7 @@
 
 #include "atlas/array/MakeView.h"
 #include "atlas/functionspace/CellColumns.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/library/config.h"
 #include "atlas/mesh/HybridElements.h"
 #include "atlas/mesh/IsGhostNode.h"
@@ -283,6 +284,8 @@ CellColumns::CellColumns(const Mesh& mesh, const eckit::Configuration& config):
         nb_cells_ = mesh.cells().size();
     }
     ATLAS_ASSERT(nb_cells_);
+
+    if (mesh_.grid()) grid_ = mesh_.grid();
 }
 
 CellColumns::~CellColumns() = default;
@@ -570,6 +573,44 @@ const parallel::Checksum& CellColumns::checksum() const {
     return *checksum_;
 }
 
+Grid CellColumns::base_grid() const {
+    if (grid_) return grid_;
+
+    const auto& comm = mpi::comm(mpi_comm());
+    std::vector<PointXY> points;
+    if (comm.size() == 1) {
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            points.push_back({xy(i, 0), xy(i, 1)});
+        }
+    } else {
+        std::vector<int> gidx;
+        std::vector<double> x, y;
+        const auto gidxView = array::make_view<gidx_t, 1>(global_index());
+        const auto ghostView = array::make_view<int, 1>(ghost());
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            if (ghostView(i) == 0) {
+                gidx.push_back(gidxView(i));
+                x.push_back(xy(i, 0));
+                y.push_back(xy(i, 1));
+            }
+        }
+        eckit::mpi::Buffer<int> gidxBuffer(comm.size());
+        eckit::mpi::Buffer<double> xBuffer(comm.size());
+        eckit::mpi::Buffer<double> yBuffer(comm.size());
+        comm.allGatherv(gidx.begin(), gidx.end(), gidxBuffer);
+        comm.allGatherv(x.begin(), x.end(), xBuffer);
+        comm.allGatherv(y.begin(), y.end(), yBuffer);
+        points.reserve(gidxBuffer.buffer.size());
+        for (auto i : gidxBuffer.buffer) {
+            points[i - 1] = atlas::PointXY{xBuffer.buffer[i - 1], yBuffer.buffer[i - 1]};
+        }
+    }
+    grid_ = UnstructuredGrid(points);
+    return grid_;
+}
+
 Field CellColumns::lonlat() const {
     if (!mesh_.cells().has_field("lonlat")) {
         mesh::actions::Build2DCellCentres("lonlat")(const_cast<Mesh&>(mesh_));
diff --git a/src/atlas/functionspace/CellColumns.h b/src/atlas/functionspace/CellColumns.h
index a53a09268..336cabba3 100644
--- a/src/atlas/functionspace/CellColumns.h
+++ b/src/atlas/functionspace/CellColumns.h
@@ -89,6 +89,8 @@ class CellColumns : public functionspace::FunctionSpaceImpl {
 
     idx_t size() const override { return nb_cells_; }
 
+    Grid base_grid() const override;
+
     Field lonlat() const override;
 
     Field remote_index() const override;
@@ -111,6 +113,7 @@ class CellColumns : public functionspace::FunctionSpaceImpl {
     size_t footprint() const override;
 
 private:                           // data
+    mutable Grid grid_;
     Mesh mesh_;                    // non-const because functionspace may modify mesh
     mesh::HybridElements& cells_;  // non-const because functionspace may modify mesh
     idx_t nb_levels_;
diff --git a/src/atlas/functionspace/EdgeColumns.cc b/src/atlas/functionspace/EdgeColumns.cc
index 07a918774..7e4e45f61 100644
--- a/src/atlas/functionspace/EdgeColumns.cc
+++ b/src/atlas/functionspace/EdgeColumns.cc
@@ -17,6 +17,7 @@
 #include "atlas/array/MakeView.h"
 #include "atlas/field/detail/FieldImpl.h"
 #include "atlas/functionspace/EdgeColumns.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/library/config.h"
 #include "atlas/mesh/HybridElements.h"
 #include "atlas/mesh/IsGhostNode.h"
@@ -281,6 +282,8 @@ EdgeColumns::EdgeColumns(const Mesh& mesh, const eckit::Configuration& config):
         nb_edges_ = get_nb_edges_from_metadata();
     }
     ATLAS_ASSERT(nb_edges_);
+
+    if (mesh_.grid()) grid_ = mesh_.grid();
 }
 
 EdgeColumns::~EdgeColumns() = default;
@@ -560,6 +563,44 @@ const parallel::Checksum& EdgeColumns::checksum() const {
     return *checksum_;
 }
 
+Grid EdgeColumns::base_grid() const {
+    if (grid_) return grid_;
+
+    const auto& comm = mpi::comm(mpi_comm());
+    std::vector<PointXY> points;
+    if (comm.size() == 1) {
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            points.push_back({xy(i, 0), xy(i, 1)});
+        }
+    } else {
+        std::vector<int> gidx;
+        std::vector<double> x, y;
+        const auto gidxView = array::make_view<gidx_t, 1>(global_index());
+        const auto ghostView = array::make_view<int, 1>(ghost());
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            if (ghostView(i) == 0) {
+                gidx.push_back(gidxView(i));
+                x.push_back(xy(i, 0));
+                y.push_back(xy(i, 1));
+            }
+        }
+        eckit::mpi::Buffer<int> gidxBuffer(comm.size());
+        eckit::mpi::Buffer<double> xBuffer(comm.size());
+        eckit::mpi::Buffer<double> yBuffer(comm.size());
+        comm.allGatherv(gidx.begin(), gidx.end(), gidxBuffer);
+        comm.allGatherv(x.begin(), x.end(), xBuffer);
+        comm.allGatherv(y.begin(), y.end(), yBuffer);
+        points.reserve(gidxBuffer.buffer.size());
+        for (auto i : gidxBuffer.buffer) {
+            points[i - 1] = atlas::PointXY{xBuffer.buffer[i - 1], yBuffer.buffer[i - 1]};
+        }
+    }
+    grid_ = UnstructuredGrid(points);
+    return grid_;
+}
+
 Field EdgeColumns::lonlat() const {
     return edges().field("lonlat");
 }
diff --git a/src/atlas/functionspace/EdgeColumns.h b/src/atlas/functionspace/EdgeColumns.h
index b3fac5a98..0140a1482 100644
--- a/src/atlas/functionspace/EdgeColumns.h
+++ b/src/atlas/functionspace/EdgeColumns.h
@@ -85,6 +85,8 @@ class EdgeColumns : public functionspace::FunctionSpaceImpl {
 
     idx_t size() const override { return nb_edges_; }
 
+    Grid base_grid() const override;
+
     Field lonlat() const override;
 
     Field global_index() const override;
@@ -105,6 +107,7 @@ class EdgeColumns : public functionspace::FunctionSpaceImpl {
     size_t footprint() const override;
 
 private:                           // data
+    mutable Grid grid_;
     Mesh mesh_;                    // non-const because functionspace may modify mesh
     mesh::HybridElements& edges_;  // non-const because functionspace may modify mesh
     idx_t nb_levels_;
diff --git a/src/atlas/functionspace/FunctionSpace.cc b/src/atlas/functionspace/FunctionSpace.cc
index cbe395c76..bdc8c60c8 100644
--- a/src/atlas/functionspace/FunctionSpace.cc
+++ b/src/atlas/functionspace/FunctionSpace.cc
@@ -11,6 +11,7 @@
 
 #include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/field/Field.h"
+#include "atlas/grid.h"
 #include "atlas/functionspace/detail/FunctionSpaceImpl.h"
 
 namespace atlas {
@@ -30,6 +31,10 @@ size_t FunctionSpace::footprint() const {
     return get()->footprint();
 }
 
+Grid FunctionSpace::base_grid() const {
+    return get()->base_grid();
+}
+
 Field FunctionSpace::createField(const eckit::Configuration& config) const {
     return get()->createField(config);
 }
diff --git a/src/atlas/functionspace/FunctionSpace.h b/src/atlas/functionspace/FunctionSpace.h
index be4e6001a..a937ef164 100644
--- a/src/atlas/functionspace/FunctionSpace.h
+++ b/src/atlas/functionspace/FunctionSpace.h
@@ -22,6 +22,7 @@ class Configuration;
 namespace atlas {
 class Field;
 class FieldSet;
+class Grid;
 class Projection;
 namespace functionspace {
 class FunctionSpaceImpl;
@@ -85,6 +86,8 @@ class FunctionSpace : DOXYGEN_HIDE(public util::ObjectHandle<functionspace::Func
 
     idx_t size() const;
 
+    Grid base_grid() const;
+
     Field lonlat() const;
 
     Field ghost() const;
diff --git a/src/atlas/functionspace/NodeColumns.cc b/src/atlas/functionspace/NodeColumns.cc
index 76fe397e3..0954f2089 100644
--- a/src/atlas/functionspace/NodeColumns.cc
+++ b/src/atlas/functionspace/NodeColumns.cc
@@ -18,6 +18,7 @@
 #include "atlas/field/FieldSet.h"
 #include "atlas/functionspace/NodeColumns.h"
 #include "atlas/grid/Grid.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/library/config.h"
 #include "atlas/mesh/IsGhostNode.h"
 #include "atlas/mesh/Mesh.h"
@@ -223,6 +224,8 @@ NodeColumns::NodeColumns(Mesh mesh, const eckit::Configuration& config):
             nb_nodes_ = mesh_.nodes().size();
         }
     }
+
+    if (mesh_.grid()) grid_ = mesh_.grid();
 }
 
 NodeColumns::~NodeColumns() = default;
@@ -641,6 +644,43 @@ const parallel::Checksum& NodeColumns::checksum() const {
 //  return checksum(fieldset);
 //}
 
+Grid NodeColumns::base_grid() const {
+    if (grid_) return grid_;
+
+    const auto& comm = mpi::comm(mpi_comm());
+    std::vector<PointXY> points;
+    if (comm.size() == 1) {
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            points.push_back({xy(i, 0), xy(i, 1)});
+        }
+    } else {
+        std::vector<int> gidx;
+        std::vector<double> x, y;
+        const auto gidxView = array::make_view<gidx_t, 1>(global_index());
+        const auto ghostView = array::make_view<int, 1>(ghost());
+        const auto xy = atlas::array::make_view<double, 2>(mesh_.nodes().xy());
+        for (auto i = 0; i < xy.shape(0); i++) {
+            if (ghostView(i) == 0) {
+                gidx.push_back(gidxView(i));
+                x.push_back(xy(i, 0));
+                y.push_back(xy(i, 1));
+            }
+        }
+        eckit::mpi::Buffer<int> gidxBuffer(comm.size());
+        eckit::mpi::Buffer<double> xBuffer(comm.size());
+        eckit::mpi::Buffer<double> yBuffer(comm.size());
+        comm.allGatherv(gidx.begin(), gidx.end(), gidxBuffer);
+        comm.allGatherv(x.begin(), x.end(), xBuffer);
+        comm.allGatherv(y.begin(), y.end(), yBuffer);
+        points.reserve(gidxBuffer.buffer.size());
+        for (auto i : gidxBuffer.buffer) {
+            points[i - 1] = atlas::PointXY{xBuffer.buffer[i - 1], yBuffer.buffer[i - 1]};
+        }
+    }
+    grid_ = UnstructuredGrid(points);
+    return grid_;
+}
 
 }  // namespace detail
 
diff --git a/src/atlas/functionspace/NodeColumns.h b/src/atlas/functionspace/NodeColumns.h
index 520d8ea71..9cc2ce989 100644
--- a/src/atlas/functionspace/NodeColumns.h
+++ b/src/atlas/functionspace/NodeColumns.h
@@ -262,6 +262,8 @@ class NodeColumns : public functionspace::FunctionSpaceImpl {
 
     idx_t nb_parts() const override { return mesh_.nb_parts(); }
 
+    Grid base_grid() const override;
+
     Field lonlat() const override { return nodes_.lonlat(); }
 
     Field ghost() const override { return nodes_.ghost(); }
@@ -293,6 +295,7 @@ class NodeColumns : public functionspace::FunctionSpaceImpl {
     virtual size_t footprint() const override { return 0; }
 
 private:                  // data
+    mutable Grid grid_;
     Mesh mesh_;           // non-const because functionspace may modify mesh
     mesh::Nodes& nodes_;  // non-const because functionspace may modify mesh
     mesh::Halo halo_;
diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc
index ffe3d542f..54df786bb 100644
--- a/src/atlas/functionspace/PointCloud.cc
+++ b/src/atlas/functionspace/PointCloud.cc
@@ -20,6 +20,7 @@
 #include "atlas/grid/Distribution.h"
 #include "atlas/grid/Grid.h"
 #include "atlas/grid/Iterator.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/option/Options.h"
 #include "atlas/parallel/HaloExchange.h"
 #include "atlas/parallel/mpi/mpi.h"
@@ -156,6 +157,7 @@ PointCloud::PointCloud(const Grid& grid, const grid::Distribution& distribution,
     auto& comm = mpi::comm(mpi_comm_);
     double halo_radius;
     config.get("halo_radius", halo_radius = 0.);
+    grid_ = grid;
 
     part_ = comm.rank();
 
@@ -278,6 +280,48 @@ Field PointCloud::ghost() const {
     return ghost_;
 }
 
+Grid PointCloud::base_grid() const {
+    if (grid_) return grid_;
+
+    std::vector<PointXY> points;
+    points.reserve(size_global_);
+    if (nb_partitions_ == 1) {
+        for (const auto& point : iterate().xy()) {
+            points.push_back(point);
+        }
+    } else {
+        std::vector<int> gidx;
+        gidx.reserve(size_global_);
+        std::vector<double> x, y;
+        x.reserve(size_global_);
+        y.reserve(size_global_);
+        const auto gidxView = array::make_view<gidx_t, 1>(global_index_);
+        const auto ghostView = array::make_view<int, 1>(ghost_);
+        int i = 0;
+        for (const auto& point : iterate().xy()) {
+            if (ghostView(i) == 0) {
+                gidx.push_back(gidxView(i));
+                x.push_back(point.x());
+                y.push_back(point.y());
+            }
+            i++;
+        }
+        const auto& comm = mpi::comm(mpi_comm_);
+        eckit::mpi::Buffer<int> gidxBuffer(comm.size());
+        eckit::mpi::Buffer<double> xBuffer(comm.size());
+        eckit::mpi::Buffer<double> yBuffer(comm.size());
+        comm.allGatherv(gidx.begin(), gidx.end(), gidxBuffer);
+        comm.allGatherv(x.begin(), x.end(), xBuffer);
+        comm.allGatherv(y.begin(), y.end(), yBuffer);
+        for (auto i : gidxBuffer.buffer) {
+            points[i - 1] = atlas::PointXY{xBuffer.buffer[i - 1], yBuffer.buffer[i - 1]};
+        }
+    }
+    grid_ = UnstructuredGrid(points);
+    return grid_;
+}
+
+
 array::ArrayShape PointCloud::config_shape(const eckit::Configuration& config) const {
     idx_t _size  = size();
     bool global(false);
diff --git a/src/atlas/functionspace/PointCloud.h b/src/atlas/functionspace/PointCloud.h
index 3e1402727..7d116085e 100644
--- a/src/atlas/functionspace/PointCloud.h
+++ b/src/atlas/functionspace/PointCloud.h
@@ -19,6 +19,7 @@
 #include "atlas/field/FieldSet.h"
 #include "atlas/functionspace/FunctionSpace.h"
 #include "atlas/functionspace/detail/FunctionSpaceImpl.h"
+#include "atlas/grid/Grid.h"
 #include "atlas/parallel/HaloExchange.h"
 #include "atlas/parallel/GatherScatter.h"
 #include "atlas/runtime/Exception.h"
@@ -57,6 +58,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl {
     operator bool() const override { return true; }
     size_t footprint() const override { return sizeof(*this); }
     std::string distribution() const override;
+    Grid base_grid() const override;
     Field lonlat() const override { return lonlat_; }
     const Field& vertical() const { return vertical_; }
     Field ghost() const override;
@@ -170,6 +172,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl {
     void create_remote_index() const;
 
 private:
+    mutable Grid grid_;
     Field lonlat_;
     Field vertical_;
     mutable Field ghost_;
diff --git a/src/atlas/functionspace/detail/BlockStructuredColumns.h b/src/atlas/functionspace/detail/BlockStructuredColumns.h
index d65ab8b59..f02878d27 100644
--- a/src/atlas/functionspace/detail/BlockStructuredColumns.h
+++ b/src/atlas/functionspace/detail/BlockStructuredColumns.h
@@ -101,6 +101,7 @@ class BlockStructuredColumns : public FunctionSpaceImpl {
     const StructuredGrid& grid() const { return structuredcolumns_->grid(); }
 
     idx_t levels() const { return structuredcolumns_->levels(); }
+    Grid base_grid() const override { return structuredcolumns_->base_grid(); }
     Field lonlat() const override { return structuredcolumns_->lonlat(); }
     Field xy() const { return structuredcolumns_->xy(); }
     Field z() const { return structuredcolumns_->z(); }
diff --git a/src/atlas/functionspace/detail/FunctionSpaceImpl.cc b/src/atlas/functionspace/detail/FunctionSpaceImpl.cc
index 36420b9fc..a882986ef 100644
--- a/src/atlas/functionspace/detail/FunctionSpaceImpl.cc
+++ b/src/atlas/functionspace/detail/FunctionSpaceImpl.cc
@@ -10,6 +10,7 @@
 
 #include "FunctionSpaceImpl.h"
 #include "atlas/field/Field.h"
+#include "atlas/grid.h"
 #include "atlas/option/Options.h"
 #include "atlas/runtime/Exception.h"
 #include "atlas/util/Metadata.h"
@@ -54,6 +55,10 @@ Field NoFunctionSpace::createField(const Field&, const eckit::Configuration&) co
     ATLAS_NOTIMPLEMENTED;
 }
 
+Grid FunctionSpaceImpl::base_grid() const {
+    ATLAS_NOTIMPLEMENTED;
+}
+
 Field FunctionSpaceImpl::lonlat() const {
     ATLAS_NOTIMPLEMENTED;
 }
diff --git a/src/atlas/functionspace/detail/FunctionSpaceImpl.h b/src/atlas/functionspace/detail/FunctionSpaceImpl.h
index f2db8d01f..f37baa28c 100644
--- a/src/atlas/functionspace/detail/FunctionSpaceImpl.h
+++ b/src/atlas/functionspace/detail/FunctionSpaceImpl.h
@@ -25,6 +25,7 @@ class Configuration;
 namespace atlas {
 class FieldSet;
 class Field;
+class Grid;
 class Projection;
 namespace util {
 class Metadata;
@@ -99,6 +100,8 @@ class FunctionSpaceImpl : public util::Object {
 
     virtual const util::PartitionPolygon& polygon(idx_t halo = 0) const;
 
+    virtual atlas::Grid base_grid() const;
+
     virtual atlas::Field lonlat() const;
 
     virtual atlas::Field ghost() const;
diff --git a/src/atlas/functionspace/detail/StructuredColumns.h b/src/atlas/functionspace/detail/StructuredColumns.h
index 07e936ea3..ce06d2c9b 100644
--- a/src/atlas/functionspace/detail/StructuredColumns.h
+++ b/src/atlas/functionspace/detail/StructuredColumns.h
@@ -117,6 +117,8 @@ class StructuredColumns : public FunctionSpaceImpl {
 
     const StructuredGrid& grid() const;
 
+    Grid base_grid() const override { return *grid_; }
+
     const Projection& projection() const override { return grid().projection(); }
 
     idx_t i_begin(idx_t j) const { return i_begin_[j]; }
diff --git a/src/tests/functionspace/test_functionspace.cc b/src/tests/functionspace/test_functionspace.cc
index 8047eb9c7..b1ebd58f1 100644
--- a/src/tests/functionspace/test_functionspace.cc
+++ b/src/tests/functionspace/test_functionspace.cc
@@ -14,11 +14,18 @@
 #include "atlas/array/MakeView.h"
 #include "atlas/field/Field.h"
 #include "atlas/field/FieldSet.h"
+#include "atlas/functionspace/BlockStructuredColumns.h"
+#include "atlas/functionspace/CellColumns.h"
+#include "atlas/functionspace/EdgeColumns.h"
 #include "atlas/functionspace/NodeColumns.h"
+#include "atlas/functionspace/PointCloud.h"
 #include "atlas/functionspace/Spectral.h"
+#include "atlas/functionspace/StructuredColumns.h"
 #include "atlas/grid/Grid.h"
 #include "atlas/grid/StructuredGrid.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/mesh/Mesh.h"
+#include "atlas/mesh/MeshBuilder.h"
 #include "atlas/mesh/Nodes.h"
 #include "atlas/meshgenerator.h"
 #include "atlas/parallel/mpi/mpi.h"
@@ -637,6 +644,73 @@ CASE("test_SpectralFunctionSpace_norm") {
         }
     }
 }
+CASE("test_functionspace_base_grid") {
+    // Create list of points and construct Grid from them
+    std::vector<PointXY> points = {
+        {0.0, 5.0}, {0.0, 0.0}, {10.0, 0.0}, {15.0, 0.0}, {5.0, 5.0}, {15.0, 5.0}
+    };
+    UnstructuredGrid grid(points);
+
+    // Create Mesh from Grid
+    Mesh mesh_from_grid(grid);
+
+    // Create Mesh from list of points using MeshBuilder
+    std::vector<double> lons;
+    std::vector<double> lats;
+    for (const auto& point : points) {
+        lons.push_back(point.x());
+        lats.push_back(point.y());
+    }
+    std::vector<int> ghosts(6, 0);
+    std::vector<gidx_t> global_indices(6);
+    std::iota(global_indices.begin(), global_indices.end(), 1);
+    const idx_t remote_index_base = 0;
+    std::vector<idx_t> remote_indices(6);
+    std::iota(remote_indices.begin(), remote_indices.end(), remote_index_base);
+    std::vector<int> partitions(6, 0);
+    std::vector<std::array<gidx_t, 3>> tri_boundary_nodes = {{{3, 6, 5}}, {{3, 4, 6}}};
+    std::vector<gidx_t> tri_global_indices                = {1, 2};
+    std::vector<std::array<gidx_t, 4>> quad_boundary_nodes = {{{1, 2, 3, 5}}};
+    std::vector<gidx_t> quad_global_indices                = {3};
+    const mesh::MeshBuilder mesh_builder{};
+    const Mesh mesh_from_meshbuilder = mesh_builder(
+        lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
+        tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices);
+
+    // Create Cell/Edge/NodeColumns and PointCloud FunctionSpaces that will save the Grid on construction
+    functionspace::CellColumns cells_from_grid(mesh_from_grid);
+    functionspace::EdgeColumns edges_from_grid(mesh_from_grid);
+    functionspace::NodeColumns nodes_from_grid(mesh_from_grid);
+    functionspace::PointCloud pointcloud_from_grid(grid);
+
+    // Create Cell/Edge/NodeColumns and PointCloud FunctionSpaces that will generate the Grid when requested
+    functionspace::CellColumns cells_from_meshbuilder(mesh_from_meshbuilder);
+    functionspace::EdgeColumns edges_from_meshbuilder(mesh_from_meshbuilder);
+    functionspace::NodeColumns nodes_from_meshbuilder(mesh_from_meshbuilder);
+    functionspace::PointCloud pointcloud_from_points(points);
+
+    // All Grids should match original
+    EXPECT(cells_from_grid.base_grid() == grid);
+    EXPECT(edges_from_grid.base_grid() == grid);
+    EXPECT(nodes_from_grid.base_grid() == grid);
+    EXPECT(pointcloud_from_grid.base_grid() == grid);
+    EXPECT(cells_from_meshbuilder.base_grid() == grid);
+    EXPECT(edges_from_meshbuilder.base_grid() == grid);
+    EXPECT(nodes_from_meshbuilder.base_grid() == grid);
+    EXPECT(pointcloud_from_points.base_grid() == grid);
+
+    // Repeat with a StructuredGrid to check StructuredColumns and BlockStructuredColumns
+    Grid structured_grid("O8");
+    Mesh structured_mesh = StructuredMeshGenerator().generate(structured_grid);
+    functionspace::StructuredColumns sc(structured_grid);
+    functionspace::BlockStructuredColumns bsc(structured_grid);
+    EXPECT(sc.base_grid() == structured_grid);
+    EXPECT(bsc.base_grid() == structured_grid);
+
+    // Check that Spectral throws exception
+    functionspace::Spectral spectral(159);
+    EXPECT_THROWS(spectral.base_grid());
+}
 
 #endif