Mesh: Remove redundant custom data pointers

For copy-on-write, we want to share attribute arrays between meshes
where possible. Mutable pointers like `Mesh.mvert` make that difficult
by making ownership vague. They also make code more complex by adding
redundancy.

The simplest solution is just removing them and retrieving layers from
`CustomData` as needed. Similar changes have already been applied to
curves and point clouds (e9f82d3dc7ee, 410a6efb747f). Removing use of
the pointers generally makes code more obvious and more reusable.

Mesh data is now accessed with a C++ API (`Mesh::edges()` or
`Mesh::edges_for_write()`), and a C API (`BKE_mesh_edges(mesh)`).

The CoW changes this commit makes possible are described in T95845
and T95842, and started in D14139 and D14140. The change also simplifies
the ongoing mesh struct-of-array refactors from T95965.

**RNA/Python Access Performance**
Theoretically, accessing mesh elements with the RNA API may become
slower, since the layer needs to be found on every random access.
However, overhead is already high enough that this doesn't make a
noticible differenc, and performance is actually improved in some
cases. Random access can be up to 10% faster, but other situations
might be a bit slower. Generally using `foreach_get/set` are the best
way to improve performance. See the differential revision for more
discussion about Python performance.

Cycles has been updated to use raw pointers and the internal Blender
mesh types, mostly because there is no sense in having this overhead
when it's already compiled with Blender. In my tests this roughly
halves the Cycles mesh creation time (0.19s to 0.10s for a 1 million
face grid).

Differential Revision: https://developer.blender.org/D15488

32 files changed, 448 insertions, 374 deletions

diff --git a/source/blender/nodes/geometry/nodes/node_geo_convex_hull.cc b/source/blender/nodes/geometry/nodes/node_geo_convex_hull.cc
index 03892501c89..d03f036eaa4 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_convex_hull.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_convex_hull.cc
@@ -46,6 +46,7 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
   }
 
   /* Copy vertices. */
+  MutableSpan<MVert> dst_verts = result->vertices_for_write();
   for (const int i : IndexRange(verts_num)) {
     float co[3];
     int original_index;
@@ -59,7 +60,7 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
       }
 #  endif
       /* Copy the position of the original point. */
-      copy_v3_v3(result->mvert[i].co, co);
+      copy_v3_v3(dst_verts[i].co, co);
     }
     else {
       BLI_assert_msg(0, "Unexpected new vertex in hull output");
@@ -73,6 +74,8 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
    * to a loop and edges need to be created from that. */
   Array<MLoop> mloop_src(loops_num);
   uint edge_index = 0;
+  MutableSpan<MEdge> edges = result->edges_for_write();
+
   for (const int i : IndexRange(loops_num)) {
     int v_from;
     int v_to;
@@ -81,7 +84,7 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
     mloop_src[i].v = (uint)v_from;
     /* Add edges for ascending order loops only. */
     if (v_from < v_to) {
-      MEdge &edge = result->medge[edge_index];
+      MEdge &edge = edges[edge_index];
       edge.v1 = v_from;
       edge.v2 = v_to;
       edge.flag = ME_EDGEDRAW | ME_EDGERENDER;
@@ -95,7 +98,7 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
   }
   if (edges_num == 1) {
     /* In this case there are no loops. */
-    MEdge &edge = result->medge[0];
+    MEdge &edge = edges[0];
     edge.v1 = 0;
     edge.v2 = 1;
     edge.flag |= ME_EDGEDRAW | ME_EDGERENDER | ME_LOOSEEDGE;
@@ -106,7 +109,10 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
   /* Copy faces. */
   Array<int> loops;
   int j = 0;
-  MLoop *loop = result->mloop;
+  MutableSpan<MPoly> polys = result->polygons_for_write();
+  MutableSpan<MLoop> mesh_loops = result->loops_for_write();
+  MLoop *loop = mesh_loops.data();
+
   for (const int i : IndexRange(faces_num)) {
     const int len = plConvexHullGetFaceSize(hull, i);
 
@@ -116,7 +122,7 @@ static Mesh *hull_from_bullet(const Mesh *mesh, Span<float3> coords)
     loops.reinitialize(len);
     plConvexHullGetFaceLoops(hull, i, loops.data());
 
-    MPoly &face = result->mpoly[i];
+    MPoly &face = polys[i];
     face.loopstart = j;
     face.totloop = len;
     for (const int k : IndexRange(len)) {
diff --git a/source/blender/nodes/geometry/nodes/node_geo_curve_fill.cc b/source/blender/nodes/geometry/nodes/node_geo_curve_fill.cc
index c9a8dba55b2..580886ad6be 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_curve_fill.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_curve_fill.cc
@@ -79,10 +79,10 @@ static Mesh *cdt_to_mesh(const meshintersect::CDT_result<double> &result)
   }
 
   Mesh *mesh = BKE_mesh_new_nomain(vert_len, edge_len, 0, loop_len, poly_len);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   for (const int i : IndexRange(result.vert.size())) {
     copy_v3_v3(verts[i].co, float3((float)result.vert[i].x, (float)result.vert[i].y, 0.0f));
diff --git a/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
index 2481170835b..07b419cc05c 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
@@ -66,6 +66,12 @@ static void deform_curves(const CurvesGeometry &curves,
 
   const float4x4 curves_to_surface = surface_to_curves.inverted();
 
+  const Span<MVert> surface_verts_old = surface_mesh_old.vertices();
+  const Span<MLoop> surface_loops_old = surface_mesh_old.loops();
+
+  const Span<MVert> surface_verts_new = surface_mesh_new.vertices();
+  const Span<MLoop> surface_loops_new = surface_mesh_new.loops();
+
   threading::parallel_for(curves.curves_range(), 256, [&](const IndexRange range) {
     for (const int curve_i : range) {
       const ReverseUVSampler::Result &surface_sample_old = surface_samples_old[curve_i];
@@ -92,13 +98,13 @@ static void deform_curves(const CurvesGeometry &curves,
       const int corner_1_new = looptri_new.tri[1];
       const int corner_2_new = looptri_new.tri[2];
 
-      const int vert_0_old = surface_mesh_old.mloop[corner_0_old].v;
-      const int vert_1_old = surface_mesh_old.mloop[corner_1_old].v;
-      const int vert_2_old = surface_mesh_old.mloop[corner_2_old].v;
+      const int vert_0_old = surface_loops_old[corner_0_old].v;
+      const int vert_1_old = surface_loops_old[corner_1_old].v;
+      const int vert_2_old = surface_loops_old[corner_2_old].v;
 
-      const int vert_0_new = surface_mesh_new.mloop[corner_0_new].v;
-      const int vert_1_new = surface_mesh_new.mloop[corner_1_new].v;
-      const int vert_2_new = surface_mesh_new.mloop[corner_2_new].v;
+      const int vert_0_new = surface_loops_new[corner_0_new].v;
+      const int vert_1_new = surface_loops_new[corner_1_new].v;
+      const int vert_2_new = surface_loops_new[corner_2_new].v;
 
       const float3 &normal_0_old = corner_normals_old[corner_0_old];
       const float3 &normal_1_old = corner_normals_old[corner_1_old];
@@ -112,14 +118,14 @@ static void deform_curves(const CurvesGeometry &curves,
       const float3 normal_new = math::normalize(
           mix3(bary_weights_new, normal_0_new, normal_1_new, normal_2_new));
 
-      const float3 &pos_0_old = surface_mesh_old.mvert[vert_0_old].co;
-      const float3 &pos_1_old = surface_mesh_old.mvert[vert_1_old].co;
-      const float3 &pos_2_old = surface_mesh_old.mvert[vert_2_old].co;
+      const float3 &pos_0_old = surface_verts_old[vert_0_old].co;
+      const float3 &pos_1_old = surface_verts_old[vert_1_old].co;
+      const float3 &pos_2_old = surface_verts_old[vert_2_old].co;
       const float3 pos_old = mix3(bary_weights_old, pos_0_old, pos_1_old, pos_2_old);
 
-      const float3 &pos_0_new = surface_mesh_new.mvert[vert_0_new].co;
-      const float3 &pos_1_new = surface_mesh_new.mvert[vert_1_new].co;
-      const float3 &pos_2_new = surface_mesh_new.mvert[vert_2_new].co;
+      const float3 &pos_0_new = surface_verts_new[vert_0_new].co;
+      const float3 &pos_1_new = surface_verts_new[vert_1_new].co;
+      const float3 &pos_2_new = surface_verts_new[vert_2_new].co;
       const float3 pos_new = mix3(bary_weights_new, pos_0_new, pos_1_new, pos_2_new);
 
       /* The translation is just the difference between the old and new position on the surface. */
diff --git a/source/blender/nodes/geometry/nodes/node_geo_delete_geometry.cc b/source/blender/nodes/geometry/nodes/node_geo_delete_geometry.cc
index 58ba2fefff9..a4096efb79f 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_delete_geometry.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_delete_geometry.cc
@@ -161,10 +161,11 @@ static void copy_face_corner_attributes(const Map<AttributeIDRef, AttributeKind>
                                         const Span<int> selected_poly_indices,
                                         const Mesh &mesh_in)
 {
+  const Span<MPoly> polys = mesh_in.polygons();
   Vector<int64_t> indices;
   indices.reserve(selected_loops_num);
   for (const int src_poly_index : selected_poly_indices) {
-    const MPoly &src_poly = mesh_in.mpoly[src_poly_index];
+    const MPoly &src_poly = polys[src_poly_index];
     const int src_loop_start = src_poly.loopstart;
     const int tot_loop = src_poly.totloop;
     for (const int i : IndexRange(tot_loop)) {
@@ -180,34 +181,30 @@ static void copy_masked_vertices_to_new_mesh(const Mesh &src_mesh,
                                              Span<int> vertex_map)
 {
   BLI_assert(src_mesh.totvert == vertex_map.size());
+  const Span<MVert> src_verts = src_mesh.vertices();
+  MutableSpan<MVert> dst_verts = dst_mesh.vertices_for_write();
+
   for (const int i_src : vertex_map.index_range()) {
     const int i_dst = vertex_map[i_src];
     if (i_dst == -1) {
       continue;
     }
-
-    const MVert &v_src = src_mesh.mvert[i_src];
-    MVert &v_dst = dst_mesh.mvert[i_dst];
-
-    v_dst = v_src;
+    dst_verts[i_dst] = src_verts[i_src];
   }
 }
 
 static void copy_masked_edges_to_new_mesh(const Mesh &src_mesh, Mesh &dst_mesh, Span<int> edge_map)
 {
   BLI_assert(src_mesh.totedge == edge_map.size());
+  const Span<MEdge> src_edges = src_mesh.edges();
+  MutableSpan<MEdge> dst_edges = dst_mesh.edges_for_write();
+
   for (const int i_src : IndexRange(src_mesh.totedge)) {
     const int i_dst = edge_map[i_src];
     if (ELEM(i_dst, -1, -2)) {
       continue;
     }
-
-    const MEdge &e_src = src_mesh.medge[i_src];
-    MEdge &e_dst = dst_mesh.medge[i_dst];
-
-    e_dst = e_src;
-    e_dst.v1 = e_src.v1;
-    e_dst.v2 = e_src.v2;
+    dst_edges[i_dst] = src_edges[i_src];
   }
 }
 
@@ -218,14 +215,16 @@ static void copy_masked_edges_to_new_mesh(const Mesh &src_mesh,
 {
   BLI_assert(src_mesh.totvert == vertex_map.size());
   BLI_assert(src_mesh.totedge == edge_map.size());
+  const Span<MEdge> src_edges = src_mesh.edges();
+  MutableSpan<MEdge> dst_edges = dst_mesh.edges_for_write();
+
   for (const int i_src : IndexRange(src_mesh.totedge)) {
     const int i_dst = edge_map[i_src];
     if (i_dst == -1) {
       continue;
     }
-
-    const MEdge &e_src = src_mesh.medge[i_src];
-    MEdge &e_dst = dst_mesh.medge[i_dst];
+    const MEdge &e_src = src_edges[i_src];
+    MEdge &e_dst = dst_edges[i_dst];
 
     e_dst = e_src;
     e_dst.v1 = vertex_map[e_src.v1];
@@ -240,16 +239,21 @@ static void copy_masked_polys_to_new_mesh(const Mesh &src_mesh,
                                           Span<int> masked_poly_indices,
                                           Span<int> new_loop_starts)
 {
+  const Span<MPoly> src_polygons = src_mesh.polygons();
+  const Span<MLoop> src_loops = src_mesh.loops();
+  MutableSpan<MPoly> dst_polygons = dst_mesh.polygons_for_write();
+  MutableSpan<MLoop> dst_loops = dst_mesh.loops_for_write();
+
   for (const int i_dst : masked_poly_indices.index_range()) {
     const int i_src = masked_poly_indices[i_dst];
 
-    const MPoly &mp_src = src_mesh.mpoly[i_src];
-    MPoly &mp_dst = dst_mesh.mpoly[i_dst];
+    const MPoly &mp_src = src_polygons[i_src];
+    MPoly &mp_dst = dst_polygons[i_dst];
     const int i_ml_src = mp_src.loopstart;
     const int i_ml_dst = new_loop_starts[i_dst];
 
-    const MLoop *ml_src = src_mesh.mloop + i_ml_src;
-    MLoop *ml_dst = dst_mesh.mloop + i_ml_dst;
+    const MLoop *ml_src = &src_loops[i_ml_src];
+    MLoop *ml_dst = &dst_loops[i_ml_dst];
 
     mp_dst = mp_src;
     mp_dst.loopstart = i_ml_dst;
@@ -266,16 +270,21 @@ static void copy_masked_polys_to_new_mesh(const Mesh &src_mesh,
                                           Span<int> masked_poly_indices,
                                           Span<int> new_loop_starts)
 {
+  const Span<MPoly> src_polygons = src_mesh.polygons();
+  const Span<MLoop> src_loops = src_mesh.loops();
+  MutableSpan<MPoly> dst_polygons = dst_mesh.polygons_for_write();
+  MutableSpan<MLoop> dst_loops = dst_mesh.loops_for_write();
+
   for (const int i_dst : masked_poly_indices.index_range()) {
     const int i_src = masked_poly_indices[i_dst];
 
-    const MPoly &mp_src = src_mesh.mpoly[i_src];
-    MPoly &mp_dst = dst_mesh.mpoly[i_dst];
+    const MPoly &mp_src = src_polygons[i_src];
+    MPoly &mp_dst = dst_polygons[i_dst];
     const int i_ml_src = mp_src.loopstart;
     const int i_ml_dst = new_loop_starts[i_dst];
 
-    const MLoop *ml_src = src_mesh.mloop + i_ml_src;
-    MLoop *ml_dst = dst_mesh.mloop + i_ml_dst;
+    const MLoop *ml_src = &src_loops[i_ml_src];
+    MLoop *ml_dst = &dst_loops[i_ml_dst];
 
     mp_dst = mp_src;
     mp_dst.loopstart = i_ml_dst;
@@ -293,16 +302,21 @@ static void copy_masked_polys_to_new_mesh(const Mesh &src_mesh,
                                           Span<int> masked_poly_indices,
                                           Span<int> new_loop_starts)
 {
+  const Span<MPoly> src_polygons = src_mesh.polygons();
+  const Span<MLoop> src_loops = src_mesh.loops();
+  MutableSpan<MPoly> dst_polygons = dst_mesh.polygons_for_write();
+  MutableSpan<MLoop> dst_loops = dst_mesh.loops_for_write();
+
   for (const int i_dst : masked_poly_indices.index_range()) {
     const int i_src = masked_poly_indices[i_dst];
 
-    const MPoly &mp_src = src_mesh.mpoly[i_src];
-    MPoly &mp_dst = dst_mesh.mpoly[i_dst];
+    const MPoly &mp_src = src_polygons[i_src];
+    MPoly &mp_dst = dst_polygons[i_dst];
     const int i_ml_src = mp_src.loopstart;
     const int i_ml_dst = new_loop_starts[i_dst];
 
-    const MLoop *ml_src = src_mesh.mloop + i_ml_src;
-    MLoop *ml_dst = dst_mesh.mloop + i_ml_dst;
+    const MLoop *ml_src = &src_loops[i_ml_src];
+    MLoop *ml_dst = &dst_loops[i_ml_dst];
 
     mp_dst = mp_src;
     mp_dst.loopstart = i_ml_dst;
@@ -419,10 +433,11 @@ static void compute_selected_edges_from_vertex_selection(const Mesh &mesh,
                                                          int *r_selected_edges_num)
 {
   BLI_assert(mesh.totedge == r_edge_map.size());
+  const Span<MEdge> edges = mesh.edges();
 
   int selected_edges_num = 0;
   for (const int i : IndexRange(mesh.totedge)) {
-    const MEdge &edge = mesh.medge[i];
+    const MEdge &edge = edges[i];
 
     /* Only add the edge if both vertices will be in the new mesh. */
     if (vertex_selection[edge.v1] && vertex_selection[edge.v2]) {
@@ -445,17 +460,19 @@ static void compute_selected_polygons_from_vertex_selection(const Mesh &mesh,
                                                             int *r_selected_loops_num)
 {
   BLI_assert(mesh.totvert == vertex_selection.size());
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
 
   r_selected_poly_indices.reserve(mesh.totpoly);
   r_loop_starts.reserve(mesh.totloop);
 
   int selected_loops_num = 0;
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly_src = mesh.mpoly[i];
+  for (const int i : polys.index_range()) {
+    const MPoly &poly_src = polys[i];
 
     bool all_verts_in_selection = true;
-    Span<MLoop> loops_src(&mesh.mloop[poly_src.loopstart], poly_src.totloop);
-    for (const MLoop &loop : loops_src) {
+    const Span<MLoop> poly_loops = loops.slice(poly_src.loopstart, poly_src.totloop);
+    for (const MLoop &loop : poly_loops) {
       if (!vertex_selection[loop.v]) {
         all_verts_in_selection = false;
         break;
@@ -486,11 +503,12 @@ static void compute_selected_vertices_and_edges_from_edge_selection(
     int *r_selected_edges_num)
 {
   BLI_assert(mesh.totedge == edge_selection.size());
+  const Span<MEdge> edges = mesh.edges();
 
   int selected_edges_num = 0;
   int selected_verts_num = 0;
   for (const int i : IndexRange(mesh.totedge)) {
-    const MEdge &edge = mesh.medge[i];
+    const MEdge &edge = edges[i];
     if (edge_selection[i]) {
       r_edge_map[i] = selected_edges_num;
       selected_edges_num++;
@@ -547,16 +565,19 @@ static void compute_selected_polygons_from_edge_selection(const Mesh &mesh,
                                                           int *r_selected_polys_num,
                                                           int *r_selected_loops_num)
 {
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   r_selected_poly_indices.reserve(mesh.totpoly);
   r_loop_starts.reserve(mesh.totloop);
 
   int selected_loops_num = 0;
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly_src = mesh.mpoly[i];
+  for (const int i : polys.index_range()) {
+    const MPoly &poly_src = polys[i];
 
     bool all_edges_in_selection = true;
-    Span<MLoop> loops_src(&mesh.mloop[poly_src.loopstart], poly_src.totloop);
-    for (const MLoop &loop : loops_src) {
+    const Span<MLoop> poly_loops = loops.slice(poly_src.loopstart, poly_src.totloop);
+    for (const MLoop &loop : poly_loops) {
       if (!edge_selection[loop.e]) {
         all_edges_in_selection = false;
         break;
@@ -654,7 +675,7 @@ static void compute_selected_mesh_data_from_edge_selection_edge_face(
 
 /**
  * Checks for every edge if it is in `edge_selection`. If it is, the vertices belonging to
- * that edge are kept as well. The polygons are kept if all edges are in the selection.
+ * that edge are kept as well. The polys are kept if all edges are in the selection.
  */
 static void compute_selected_mesh_data_from_edge_selection(const Mesh &mesh,
                                                            const Span<bool> edge_selection,
@@ -693,13 +714,14 @@ static void compute_selected_polygons_from_poly_selection(const Mesh &mesh,
                                                           int *r_selected_loops_num)
 {
   BLI_assert(mesh.totpoly == poly_selection.size());
+  const Span<MPoly> polys = mesh.polygons();
 
   r_selected_poly_indices.reserve(mesh.totpoly);
   r_loop_starts.reserve(mesh.totloop);
 
   int selected_loops_num = 0;
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly_src = mesh.mpoly[i];
+  for (const int i : polys.index_range()) {
+    const MPoly &poly_src = polys[i];
     /* We keep this one. */
     if (poly_selection[i]) {
       r_selected_poly_indices.append_unchecked(i);
@@ -726,6 +748,9 @@ static void compute_selected_mesh_data_from_poly_selection_edge_face(
 {
   BLI_assert(mesh.totpoly == poly_selection.size());
   BLI_assert(mesh.totedge == r_edge_map.size());
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   r_edge_map.fill(-1);
 
   r_selected_poly_indices.reserve(mesh.totpoly);
@@ -733,8 +758,8 @@ static void compute_selected_mesh_data_from_poly_selection_edge_face(
 
   int selected_loops_num = 0;
   int selected_edges_num = 0;
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly_src = mesh.mpoly[i];
+  for (const int i : polys.index_range()) {
+    const MPoly &poly_src = polys[i];
     /* We keep this one. */
     if (poly_selection[i]) {
       r_selected_poly_indices.append_unchecked(i);
@@ -742,8 +767,8 @@ static void compute_selected_mesh_data_from_poly_selection_edge_face(
       selected_loops_num += poly_src.totloop;
 
       /* Add the vertices and the edges. */
-      Span<MLoop> loops_src(&mesh.mloop[poly_src.loopstart], poly_src.totloop);
-      for (const MLoop &loop : loops_src) {
+      const Span<MLoop> poly_loops = loops.slice(poly_src.loopstart, poly_src.totloop);
+      for (const MLoop &loop : poly_loops) {
         /* Check first if it has not yet been added. */
         if (r_edge_map[loop.e] == -1) {
           r_edge_map[loop.e] = selected_edges_num;
@@ -774,6 +799,9 @@ static void compute_selected_mesh_data_from_poly_selection(const Mesh &mesh,
 {
   BLI_assert(mesh.totpoly == poly_selection.size());
   BLI_assert(mesh.totedge == r_edge_map.size());
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   r_vertex_map.fill(-1);
   r_edge_map.fill(-1);
 
@@ -783,8 +811,8 @@ static void compute_selected_mesh_data_from_poly_selection(const Mesh &mesh,
   int selected_loops_num = 0;
   int selected_verts_num = 0;
   int selected_edges_num = 0;
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly_src = mesh.mpoly[i];
+  for (const int i : polys.index_range()) {
+    const MPoly &poly_src = polys[i];
     /* We keep this one. */
     if (poly_selection[i]) {
       r_selected_poly_indices.append_unchecked(i);
@@ -792,8 +820,8 @@ static void compute_selected_mesh_data_from_poly_selection(const Mesh &mesh,
       selected_loops_num += poly_src.totloop;
 
       /* Add the vertices and the edges. */
-      Span<MLoop> loops_src(&mesh.mloop[poly_src.loopstart], poly_src.totloop);
-      for (const MLoop &loop : loops_src) {
+      const Span<MLoop> poly_loops = loops.slice(poly_src.loopstart, poly_src.totloop);
+      for (const MLoop &loop : poly_loops) {
         /* Check first if it has not yet been added. */
         if (r_vertex_map[loop.v] == -1) {
           r_vertex_map[loop.v] = selected_verts_num;
@@ -968,7 +996,7 @@ static void do_mesh_separation(GeometrySet &geometry_set,
                                                    selected_polys_num);
 
       /* Copy the selected parts of the mesh over to the new mesh. */
-      memcpy(mesh_out->mvert, mesh_in.mvert, mesh_in.totvert * sizeof(MVert));
+      mesh_out->vertices_for_write().copy_from(mesh_in.vertices());
       copy_masked_edges_to_new_mesh(mesh_in, *mesh_out, edge_map);
       copy_masked_polys_to_new_mesh(
           mesh_in, *mesh_out, edge_map, selected_poly_indices, new_loop_starts);
@@ -1031,8 +1059,8 @@ static void do_mesh_separation(GeometrySet &geometry_set,
           &mesh_in, mesh_in.totvert, mesh_in.totedge, 0, selected_loops_num, selected_polys_num);
 
       /* Copy the selected parts of the mesh over to the new mesh. */
-      memcpy(mesh_out->mvert, mesh_in.mvert, mesh_in.totvert * sizeof(MVert));
-      memcpy(mesh_out->medge, mesh_in.medge, mesh_in.totedge * sizeof(MEdge));
+      mesh_out->vertices_for_write().copy_from(mesh_in.vertices());
+      mesh_out->edges_for_write().copy_from(mesh_in.edges());
       copy_masked_polys_to_new_mesh(mesh_in, *mesh_out, selected_poly_indices, new_loop_starts);
 
       /* Copy attributes. */
diff --git a/source/blender/nodes/geometry/nodes/node_geo_distribute_points_on_faces.cc b/source/blender/nodes/geometry/nodes/node_geo_distribute_points_on_faces.cc
index d9115d39705..2bb321a349a 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_distribute_points_on_faces.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_distribute_points_on_faces.cc
@@ -105,6 +105,8 @@ static void sample_mesh_surface(const Mesh &mesh,
                                 Vector<float3> &r_bary_coords,
                                 Vector<int> &r_looptri_indices)
 {
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MLoop> loops = mesh.loops();
   const Span<MLoopTri> looptris{BKE_mesh_runtime_looptri_ensure(&mesh),
                                 BKE_mesh_runtime_looptri_len(&mesh)};
 
@@ -113,12 +115,12 @@ static void sample_mesh_surface(const Mesh &mesh,
     const int v0_loop = looptri.tri[0];
     const int v1_loop = looptri.tri[1];
     const int v2_loop = looptri.tri[2];
-    const int v0_index = mesh.mloop[v0_loop].v;
-    const int v1_index = mesh.mloop[v1_loop].v;
-    const int v2_index = mesh.mloop[v2_loop].v;
-    const float3 v0_pos = float3(mesh.mvert[v0_index].co);
-    const float3 v1_pos = float3(mesh.mvert[v1_index].co);
-    const float3 v2_pos = float3(mesh.mvert[v2_index].co);
+    const int v0_index = loops[v0_loop].v;
+    const int v1_index = loops[v1_loop].v;
+    const int v2_index = loops[v2_loop].v;
+    const float3 v0_pos = verts[v0_index].co;
+    const float3 v1_pos = verts[v1_index].co;
+    const float3 v2_pos = verts[v2_index].co;
 
     float looptri_density_factor = 1.0f;
     if (!density_factors.is_empty()) {
@@ -348,6 +350,8 @@ BLI_NOINLINE static void compute_attribute_outputs(const Mesh &mesh,
         attribute_outputs.rotation_id.get(), ATTR_DOMAIN_POINT);
   }
 
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MLoop> loops = mesh.loops();
   const Span<MLoopTri> looptris{BKE_mesh_runtime_looptri_ensure(&mesh),
                                 BKE_mesh_runtime_looptri_len(&mesh)};
 
@@ -356,12 +360,12 @@ BLI_NOINLINE static void compute_attribute_outputs(const Mesh &mesh,
     const MLoopTri &looptri = looptris[looptri_index];
     const float3 &bary_coord = bary_coords[i];
 
-    const int v0_index = mesh.mloop[looptri.tri[0]].v;
-    const int v1_index = mesh.mloop[looptri.tri[1]].v;
-    const int v2_index = mesh.mloop[looptri.tri[2]].v;
-    const float3 v0_pos = float3(mesh.mvert[v0_index].co);
-    const float3 v1_pos = float3(mesh.mvert[v1_index].co);
-    const float3 v2_pos = float3(mesh.mvert[v2_index].co);
+    const int v0_index = loops[looptri.tri[0]].v;
+    const int v1_index = loops[looptri.tri[1]].v;
+    const int v2_index = loops[looptri.tri[2]].v;
+    const float3 v0_pos = verts[v0_index].co;
+    const float3 v1_pos = verts[v1_index].co;
+    const float3 v2_pos = verts[v2_index].co;
 
     ids.span[i] = noise::hash(noise::hash_float(bary_coord), looptri_index);
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_dual_mesh.cc b/source/blender/nodes/geometry/nodes/node_geo_dual_mesh.cc
index 76eeee95239..1b9e9ae9b4a 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_dual_mesh.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_dual_mesh.cc
@@ -209,13 +209,18 @@ static void calc_boundaries(const Mesh &mesh,
 {
   BLI_assert(r_vertex_types.size() == mesh.totvert);
   BLI_assert(r_edge_types.size() == mesh.totedge);
+  const Span<MEdge> edges = mesh.edges();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   r_vertex_types.fill(VertexType::Loose);
   r_edge_types.fill(EdgeType::Loose);
 
   /* Add up the number of polys connected to each edge. */
   for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly = mesh.mpoly[i];
-    for (const MLoop &loop : Span<MLoop>(&mesh.mloop[poly.loopstart], poly.totloop)) {
+    const MPoly &poly = polys[i];
+    const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+    for (const MLoop &loop : poly_loops) {
       r_edge_types[loop.e] = get_edge_type_with_added_neighbor(r_edge_types[loop.e]);
     }
   }
@@ -226,7 +231,7 @@ static void calc_boundaries(const Mesh &mesh,
     if (edge_type == EdgeType::Loose) {
       continue;
     }
-    const MEdge &edge = mesh.medge[i];
+    const MEdge &edge = edges[i];
     if (edge_type == EdgeType::Boundary) {
       r_vertex_types[edge.v1] = get_vertex_type_with_added_neighbor(r_vertex_types[edge.v1]);
       r_vertex_types[edge.v2] = get_vertex_type_with_added_neighbor(r_vertex_types[edge.v2]);
@@ -241,7 +246,7 @@ static void calc_boundaries(const Mesh &mesh,
   for (const int i : IndexRange(mesh.totedge)) {
     const EdgeType edge_type = r_edge_types[i];
     if (edge_type == EdgeType::Normal) {
-      const MEdge &edge = mesh.medge[i];
+      const MEdge &edge = edges[i];
       if (r_vertex_types[edge.v1] == VertexType::Loose) {
         r_vertex_types[edge.v1] = VertexType::Normal;
       }
@@ -258,9 +263,12 @@ static void calc_boundaries(const Mesh &mesh,
 static void create_vertex_poly_map(const Mesh &mesh,
                                    MutableSpan<Vector<int>> r_vertex_poly_indices)
 {
-  for (const int i : IndexRange(mesh.totpoly)) {
-    const MPoly &poly = mesh.mpoly[i];
-    for (const MLoop &loop : Span<MLoop>(&mesh.mloop[poly.loopstart], poly.totloop)) {
+  const Span<MPoly> polygons = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+  for (const int i : polygons.index_range()) {
+    const MPoly &poly = polygons[i];
+    const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+    for (const MLoop &loop : poly_loops) {
       r_vertex_poly_indices[loop.v].append(i);
     }
   }
@@ -321,7 +329,9 @@ static void create_vertex_poly_map(const Mesh &mesh,
  * - Finally if we are in the normal case we also need to add the last "shared edge" to close the
  *   loop.
  */
-static bool sort_vertex_polys(const Mesh &mesh,
+static bool sort_vertex_polys(const Span<MEdge> edges,
+                              const Span<MPoly> polys,
+                              const Span<MLoop> loops,
                               const int vertex_index,
                               const bool boundary_vertex,
                               const Span<EdgeType> edge_types,
@@ -336,11 +346,11 @@ static bool sort_vertex_polys(const Mesh &mesh,
   /* For each polygon store the two corners whose edge contains the vertex. */
   Array<std::pair<int, int>> poly_vertex_corners(connected_polygons.size());
   for (const int i : connected_polygons.index_range()) {
-    const MPoly &poly = mesh.mpoly[connected_polygons[i]];
+    const MPoly &poly = polys[connected_polygons[i]];
     bool first_edge_done = false;
     for (const int loop_index : IndexRange(poly.loopstart, poly.totloop)) {
-      const MLoop &loop = mesh.mloop[loop_index];
-      if (mesh.medge[loop.e].v1 == vertex_index || mesh.medge[loop.e].v2 == vertex_index) {
+      const MLoop &loop = loops[loop_index];
+      if (edges[loop.e].v1 == vertex_index || edges[loop.e].v2 == vertex_index) {
         if (!first_edge_done) {
           poly_vertex_corners[i].first = loop_index;
           first_edge_done = true;
@@ -360,8 +370,8 @@ static bool sort_vertex_polys(const Mesh &mesh,
   if (boundary_vertex) {
     /* Our first polygon needs to be one which has a boundary edge. */
     for (const int i : connected_polygons.index_range()) {
-      const MLoop &first_loop = mesh.mloop[poly_vertex_corners[i].first];
-      const MLoop &second_loop = mesh.mloop[poly_vertex_corners[i].second];
+      const MLoop &first_loop = loops[poly_vertex_corners[i].first];
+      const MLoop &second_loop = loops[poly_vertex_corners[i].second];
       if (edge_types[first_loop.e] == EdgeType::Boundary && first_loop.v == vertex_index) {
         shared_edge_i = second_loop.e;
         r_sorted_corners[0] = poly_vertex_corners[i].first;
@@ -381,8 +391,8 @@ static bool sort_vertex_polys(const Mesh &mesh,
       /* The rotation is inconsistent between the two polygons on the boundary. Just choose one
        * of the polygon's orientation. */
       for (const int i : connected_polygons.index_range()) {
-        const MLoop &first_loop = mesh.mloop[poly_vertex_corners[i].first];
-        const MLoop &second_loop = mesh.mloop[poly_vertex_corners[i].second];
+        const MLoop &first_loop = loops[poly_vertex_corners[i].first];
+        const MLoop &second_loop = loops[poly_vertex_corners[i].second];
         if (edge_types[first_loop.e] == EdgeType::Boundary) {
           shared_edge_i = second_loop.e;
           r_sorted_corners[0] = poly_vertex_corners[i].first;
@@ -402,8 +412,8 @@ static bool sort_vertex_polys(const Mesh &mesh,
   }
   else {
     /* Any polygon can be the first. Just need to check the orientation. */
-    const MLoop &first_loop = mesh.mloop[poly_vertex_corners[0].first];
-    const MLoop &second_loop = mesh.mloop[poly_vertex_corners[0].second];
+    const MLoop &first_loop = loops[poly_vertex_corners[0].first];
+    const MLoop &second_loop = loops[poly_vertex_corners[0].second];
     if (first_loop.v == vertex_index) {
       shared_edge_i = second_loop.e;
       r_sorted_corners[0] = poly_vertex_corners[0].first;
@@ -421,8 +431,8 @@ static bool sort_vertex_polys(const Mesh &mesh,
     /* Look at the other polys to see if it has this shared edge. */
     int j = i + 1;
     for (; j < connected_polygons.size(); ++j) {
-      const MLoop &first_loop = mesh.mloop[poly_vertex_corners[j].first];
-      const MLoop &second_loop = mesh.mloop[poly_vertex_corners[j].second];
+      const MLoop &first_loop = loops[poly_vertex_corners[j].first];
+      const MLoop &second_loop = loops[poly_vertex_corners[j].second];
       if (first_loop.e == shared_edge_i) {
         r_sorted_corners[i + 1] = poly_vertex_corners[j].first;
         shared_edge_i = second_loop.e;
@@ -455,14 +465,16 @@ static bool sort_vertex_polys(const Mesh &mesh,
  * Get the edge on the poly that contains the given vertex and is a boundary edge.
  */
 static void boundary_edge_on_poly(const MPoly &poly,
-                                  const Mesh &mesh,
+                                  const Span<MEdge> edges,
+                                  const Span<MLoop> loops,
                                   const int vertex_index,
                                   const Span<EdgeType> edge_types,
                                   int &r_edge)
 {
-  for (const MLoop &loop : Span<MLoop>(&mesh.mloop[poly.loopstart], poly.totloop)) {
+  const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+  for (const MLoop &loop : poly_loops) {
     if (edge_types[loop.e] == EdgeType::Boundary) {
-      const MEdge &edge = mesh.medge[loop.e];
+      const MEdge &edge = edges[loop.e];
       if (edge.v1 == vertex_index || edge.v2 == vertex_index) {
         r_edge = loop.e;
         return;
@@ -476,7 +488,8 @@ static void boundary_edge_on_poly(const MPoly &poly,
  * orientation of the poly is taken into account.
  */
 static void boundary_edges_on_poly(const MPoly &poly,
-                                   const Mesh &mesh,
+                                   const Span<MEdge> edges,
+                                   const Span<MLoop> loops,
                                    const int vertex_index,
                                    const Span<EdgeType> edge_types,
                                    int &r_edge1,
@@ -486,9 +499,10 @@ static void boundary_edges_on_poly(const MPoly &poly,
   /* This is set to true if the order in which we encounter the two edges is inconsistent with the
    * orientation of the polygon. */
   bool needs_swap = false;
-  for (const MLoop &loop : Span<MLoop>(&mesh.mloop[poly.loopstart], poly.totloop)) {
+  const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+  for (const MLoop &loop : poly_loops) {
     if (edge_types[loop.e] == EdgeType::Boundary) {
-      const MEdge &edge = mesh.medge[loop.e];
+      const MEdge &edge = edges[loop.e];
       if (edge.v1 == vertex_index || edge.v2 == vertex_index) {
         if (edge1_done) {
           if (needs_swap) {
@@ -510,7 +524,7 @@ static void boundary_edges_on_poly(const MPoly &poly,
   }
 }
 
-static void add_edge(const Mesh &mesh,
+static void add_edge(const Span<MEdge> src_edges,
                      const int old_edge_i,
                      const int v1,
                      const int v2,
@@ -518,7 +532,7 @@ static void add_edge(const Mesh &mesh,
                      Vector<MEdge> &new_edges,
                      Vector<int> &loop_edges)
 {
-  MEdge new_edge = MEdge(mesh.medge[old_edge_i]);
+  MEdge new_edge = src_edges[old_edge_i];
   new_edge.v1 = v1;
   new_edge.v2 = v2;
   const int new_edge_i = new_edges.size();
@@ -549,14 +563,17 @@ static bool vertex_needs_dissolving(const int vertex,
  * edges in the input mesh which contain such a vertex are marked as 'done' to prevent duplicate
  * edges being created. (See T94144)
  */
-static void dissolve_redundant_verts(const Mesh &mesh,
+static void dissolve_redundant_verts(const Span<MEdge> edges,
+                                     const Span<MPoly> polys,
+                                     const Span<MLoop> loops,
                                      const Span<Vector<int>> vertex_poly_indices,
                                      MutableSpan<VertexType> vertex_types,
                                      MutableSpan<int> old_to_new_edges_map,
                                      Vector<MEdge> &new_edges,
                                      Vector<int> &new_to_old_edges_map)
 {
-  for (const int vert_i : IndexRange(mesh.totvert)) {
+  const int vertex_num = vertex_types.size();
+  for (const int vert_i : IndexRange(vertex_num)) {
     if (vertex_poly_indices[vert_i].size() != 2 || vertex_types[vert_i] != VertexType::Normal) {
       continue;
     }
@@ -564,9 +581,10 @@ static void dissolve_redundant_verts(const Mesh &mesh,
     const int second_poly_index = vertex_poly_indices[vert_i][1];
     const int new_edge_index = new_edges.size();
     bool edge_created = false;
-    const MPoly &poly = mesh.mpoly[first_poly_index];
-    for (const MLoop &loop : Span<MLoop>(&mesh.mloop[poly.loopstart], poly.totloop)) {
-      const MEdge &edge = mesh.medge[loop.e];
+    const MPoly &poly = polys[first_poly_index];
+    const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
+    for (const MLoop &loop : poly_loops) {
+      const MEdge &edge = edges[loop.e];
       const int v1 = edge.v1;
       const int v2 = edge.v2;
       bool mark_edge = false;
@@ -617,6 +635,10 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
                            const bool keep_boundaries)
 {
   const Mesh &mesh_in = *in_component.get_for_read();
+  const Span<MVert> src_verts = mesh_in.vertices();
+  const Span<MEdge> src_edges = mesh_in.edges();
+  const Span<MPoly> src_polys = mesh_in.polygons();
+  const Span<MLoop> src_loops = mesh_in.loops();
 
   Map<AttributeIDRef, AttributeKind> attributes;
   geometry_set.gather_attributes_for_propagation(
@@ -644,14 +666,28 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
       bool vertex_ok = true;
       if (vertex_types[i] == VertexType::Normal) {
         Array<int> shared_edges(loop_indices.size());
-        vertex_ok = sort_vertex_polys(
-            mesh_in, i, false, edge_types, loop_indices, shared_edges, sorted_corners);
+        vertex_ok = sort_vertex_polys(src_edges,
+                                      src_polys,
+                                      src_loops,
+                                      i,
+                                      false,
+                                      edge_types,
+                                      loop_indices,
+                                      shared_edges,
+                                      sorted_corners);
         vertex_shared_edges[i] = std::move(shared_edges);
       }
       else {
         Array<int> shared_edges(loop_indices.size() - 1);
-        vertex_ok = sort_vertex_polys(
-            mesh_in, i, true, edge_types, loop_indices, shared_edges, sorted_corners);
+        vertex_ok = sort_vertex_polys(src_edges,
+                                      src_polys,
+                                      src_loops,
+                                      i,
+                                      true,
+                                      edge_types,
+                                      loop_indices,
+                                      shared_edges,
+                                      sorted_corners);
         vertex_shared_edges[i] = std::move(shared_edges);
       }
       if (!vertex_ok) {
@@ -666,9 +702,9 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
 
   Vector<float3> vertex_positions(mesh_in.totpoly);
   for (const int i : IndexRange(mesh_in.totpoly)) {
-    const MPoly poly = mesh_in.mpoly[i];
+    const MPoly &poly = src_polys[i];
     BKE_mesh_calc_poly_center(
-        &poly, &mesh_in.mloop[poly.loopstart], mesh_in.mvert, vertex_positions[i]);
+        &poly, &src_loops[poly.loopstart], src_verts.data(), vertex_positions[i]);
   }
 
   Array<int> boundary_edge_midpoint_index;
@@ -679,8 +715,8 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
     for (const int i : IndexRange(mesh_in.totedge)) {
       if (edge_types[i] == EdgeType::Boundary) {
         float3 mid;
-        const MEdge &edge = mesh_in.medge[i];
-        mid_v3_v3v3(mid, mesh_in.mvert[edge.v1].co, mesh_in.mvert[edge.v2].co);
+        const MEdge &edge = src_edges[i];
+        mid_v3_v3v3(mid, src_verts[edge.v1].co, src_verts[edge.v2].co);
         boundary_edge_midpoint_index[i] = vertex_positions.size();
         vertex_positions.append(mid);
       }
@@ -706,7 +742,9 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
 
   /* This is necessary to prevent duplicate edges from being created, but will likely not do
    * anything for most meshes. */
-  dissolve_redundant_verts(mesh_in,
+  dissolve_redundant_verts(src_edges,
+                           src_polys,
+                           src_loops,
                            vertex_poly_indices,
                            vertex_types,
                            old_to_new_edges_map,
@@ -734,7 +772,7 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
         const int old_edge_i = shared_edges[i];
         if (old_to_new_edges_map[old_edge_i] == -1) {
           /* This edge has not been created yet. */
-          MEdge new_edge = MEdge(mesh_in.medge[old_edge_i]);
+          MEdge new_edge = src_edges[old_edge_i];
           new_edge.v1 = loop_indices[i];
           new_edge.v2 = loop_indices[(i + 1) % loop_indices.size()];
           new_to_old_edges_map.append(old_edge_i);
@@ -776,7 +814,7 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
         const int old_edge_i = shared_edges[i];
         if (old_to_new_edges_map[old_edge_i] == -1) {
           /* This edge has not been created yet. */
-          MEdge new_edge = MEdge(mesh_in.medge[old_edge_i]);
+          MEdge new_edge = src_edges[old_edge_i];
           new_edge.v1 = loop_indices[i];
           new_edge.v2 = loop_indices[i + 1];
           new_to_old_edges_map.append(old_edge_i);
@@ -795,13 +833,15 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
       int edge2;
       if (loop_indices.size() >= 2) {
         /* The first boundary edge is at the end of the chain of polygons. */
-        boundary_edge_on_poly(mesh_in.mpoly[loop_indices.last()], mesh_in, i, edge_types, edge1);
-        boundary_edge_on_poly(mesh_in.mpoly[loop_indices.first()], mesh_in, i, edge_types, edge2);
+        boundary_edge_on_poly(
+            src_polys[loop_indices.last()], src_edges, src_loops, i, edge_types, edge1);
+        boundary_edge_on_poly(
+            src_polys[loop_indices.first()], src_edges, src_loops, i, edge_types, edge2);
       }
       else {
         /* If there is only one polygon both edges are in that polygon. */
         boundary_edges_on_poly(
-            mesh_in.mpoly[loop_indices[0]], mesh_in, i, edge_types, edge1, edge2);
+            src_polys[loop_indices[0]], src_edges, src_loops, i, edge_types, edge1, edge2);
       }
 
       const int last_face_center = loop_indices.last();
@@ -809,7 +849,7 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
       new_to_old_face_corners_map.append(sorted_corners.last());
       const int first_midpoint = loop_indices.last();
       if (old_to_new_edges_map[edge1] == -1) {
-        add_edge(mesh_in,
+        add_edge(src_edges,
                  edge1,
                  last_face_center,
                  first_midpoint,
@@ -827,9 +867,9 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
       new_to_old_face_corners_map.append(sorted_corners.first());
       boundary_vertex_to_relevant_face_map.append(
           std::pair(loop_indices.last(), last_face_center));
-      vertex_positions.append(mesh_in.mvert[i].co);
+      vertex_positions.append(src_verts[i].co);
       const int boundary_vertex = loop_indices.last();
-      add_edge(mesh_in,
+      add_edge(src_edges,
                edge1,
                first_midpoint,
                boundary_vertex,
@@ -840,7 +880,7 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
       loop_indices.append(boundary_edge_midpoint_index[edge2]);
       new_to_old_face_corners_map.append(sorted_corners.first());
       const int second_midpoint = loop_indices.last();
-      add_edge(mesh_in,
+      add_edge(src_edges,
                edge2,
                boundary_vertex,
                second_midpoint,
@@ -850,7 +890,7 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
 
       if (old_to_new_edges_map[edge2] == -1) {
         const int first_face_center = loop_indices.first();
-        add_edge(mesh_in,
+        add_edge(src_edges,
                  edge2,
                  second_midpoint,
                  first_face_center,
@@ -881,20 +921,25 @@ static void calc_dual_mesh(GeometrySet &geometry_set,
                       bke::mesh_attributes(mesh_in),
                       bke::mesh_attributes_for_write(*mesh_out));
 
+  MutableSpan<MVert> dst_verts = mesh_out->vertices_for_write();
+  MutableSpan<MEdge> dst_edges = mesh_out->edges_for_write();
+  MutableSpan<MPoly> dst_polys = mesh_out->polygons_for_write();
+  MutableSpan<MLoop> dst_loops = mesh_out->loops_for_write();
+
   int loop_start = 0;
   for (const int i : IndexRange(mesh_out->totpoly)) {
-    mesh_out->mpoly[i].loopstart = loop_start;
-    mesh_out->mpoly[i].totloop = loop_lengths[i];
+    dst_polys[i].loopstart = loop_start;
+    dst_polys[i].totloop = loop_lengths[i];
     loop_start += loop_lengths[i];
   }
   for (const int i : IndexRange(mesh_out->totloop)) {
-    mesh_out->mloop[i].v = loops[i];
-    mesh_out->mloop[i].e = loop_edges[i];
+    dst_loops[i].v = loops[i];
+    dst_loops[i].e = loop_edges[i];
   }
   for (const int i : IndexRange(mesh_out->totvert)) {
-    copy_v3_v3(mesh_out->mvert[i].co, vertex_positions[i]);
+    copy_v3_v3(dst_verts[i].co, vertex_positions[i]);
   }
-  memcpy(mesh_out->medge, new_edges.data(), sizeof(MEdge) * new_edges.size());
+  dst_edges.copy_from(new_edges);
   geometry_set.replace_mesh(mesh_out);
 }
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_duplicate_elements.cc b/source/blender/nodes/geometry/nodes/node_geo_duplicate_elements.cc
index 2eb3706bac9..9af1536a194 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_duplicate_elements.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_duplicate_elements.cc
@@ -486,7 +486,7 @@ static void copy_stable_id_faces(const Mesh &mesh,
   VArraySpan<int> src{src_attribute.varray.typed<int>()};
   MutableSpan<int> dst = dst_attribute.span.typed<int>();
 
-  Span<MPoly> polys(mesh.mpoly, mesh.totpoly);
+  const Span<MPoly> polys = mesh.polygons();
   int loop_index = 0;
   for (const int i_poly : selection.index_range()) {
     const IndexRange range = range_for_offsets_index(poly_offsets, i_poly);
@@ -522,10 +522,10 @@ static void duplicate_faces(GeometrySet &geometry_set,
   geometry_set.keep_only_during_modify({GEO_COMPONENT_TYPE_MESH});
 
   const Mesh &mesh = *geometry_set.get_mesh_for_read();
-  Span<MVert> verts(mesh.mvert, mesh.totvert);
-  Span<MEdge> edges(mesh.medge, mesh.totedge);
-  Span<MPoly> polys(mesh.mpoly, mesh.totpoly);
-  Span<MLoop> loops(mesh.mloop, mesh.totloop);
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MEdge> edges = mesh.edges();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
 
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_FACE};
   FieldEvaluator evaluator(field_context, polys.size());
@@ -547,10 +547,10 @@ static void duplicate_faces(GeometrySet &geometry_set,
   offsets[selection.size()] = total_polys;
 
   Mesh *new_mesh = BKE_mesh_new_nomain(total_loops, total_loops, 0, total_loops, total_polys);
-  MutableSpan<MVert> new_verts(new_mesh->mvert, new_mesh->totvert);
-  MutableSpan<MEdge> new_edges(new_mesh->medge, new_mesh->totedge);
-  MutableSpan<MLoop> new_loops(new_mesh->mloop, new_mesh->totloop);
-  MutableSpan<MPoly> new_poly(new_mesh->mpoly, new_mesh->totpoly);
+  MutableSpan<MVert> new_verts = new_mesh->vertices_for_write();
+  MutableSpan<MEdge> new_edges = new_mesh->edges_for_write();
+  MutableSpan<MPoly> new_polys = new_mesh->polygons_for_write();
+  MutableSpan<MLoop> new_loops = new_mesh->loops_for_write();
 
   Array<int> vert_mapping(new_verts.size());
   Array<int> edge_mapping(new_edges.size());
@@ -563,8 +563,8 @@ static void duplicate_faces(GeometrySet &geometry_set,
 
     const MPoly &source = polys[selection[i_selection]];
     for ([[maybe_unused]] const int i_duplicate : IndexRange(poly_range.size())) {
-      new_poly[poly_index] = source;
-      new_poly[poly_index].loopstart = loop_index;
+      new_polys[poly_index] = source;
+      new_polys[poly_index].loopstart = loop_index;
       for (const int i_loops : IndexRange(source.totloop)) {
         const MLoop &current_loop = loops[source.loopstart + i_loops];
         loop_mapping[loop_index] = source.loopstart + i_loops;
@@ -577,7 +577,7 @@ static void duplicate_faces(GeometrySet &geometry_set,
           new_edges[loop_index].v2 = loop_index + 1;
         }
         else {
-          new_edges[loop_index].v2 = new_poly[poly_index].loopstart;
+          new_edges[loop_index].v2 = new_polys[poly_index].loopstart;
         }
         new_loops[loop_index].v = loop_index;
         new_loops[loop_index].e = loop_index;
@@ -689,7 +689,7 @@ static void copy_stable_id_edges(const Mesh &mesh,
     return;
   }
 
-  Span<MEdge> edges(mesh.medge, mesh.totedge);
+  const Span<MEdge> edges = mesh.edges();
 
   VArraySpan<int> src{src_attribute.varray.typed<int>()};
   MutableSpan<int> dst = dst_attribute.span.typed<int>();
@@ -723,8 +723,7 @@ static void duplicate_edges(GeometrySet &geometry_set,
     return;
   };
   const Mesh &mesh = *geometry_set.get_mesh_for_read();
-  Span<MVert> verts(mesh.mvert, mesh.totvert);
-  Span<MEdge> edges(mesh.medge, mesh.totedge);
+  const Span<MEdge> edges = mesh.edges();
 
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_EDGE};
   FieldEvaluator evaluator{field_context, edges.size()};
@@ -737,8 +736,7 @@ static void duplicate_edges(GeometrySet &geometry_set,
   Array<int> edge_offsets = accumulate_counts_to_offsets(selection, counts);
 
   Mesh *new_mesh = BKE_mesh_new_nomain(edge_offsets.last() * 2, edge_offsets.last(), 0, 0, 0);
-  MutableSpan<MVert> new_verts(new_mesh->mvert, new_mesh->totvert);
-  MutableSpan<MEdge> new_edges(new_mesh->medge, new_mesh->totedge);
+  MutableSpan<MEdge> new_edges = new_mesh->edges_for_write();
 
   Array<int> vert_orig_indices(edge_offsets.last() * 2);
   threading::parallel_for(selection.index_range(), 1024, [&](IndexRange range) {
@@ -906,7 +904,7 @@ static void duplicate_points_mesh(GeometrySet &geometry_set,
                                   const IndexAttributes &attribute_outputs)
 {
   const Mesh &mesh = *geometry_set.get_mesh_for_read();
-  Span<MVert> src_verts(mesh.mvert, mesh.totvert);
+  const Span<MVert> src_verts = mesh.vertices();
 
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_POINT};
   FieldEvaluator evaluator{field_context, src_verts.size()};
@@ -919,7 +917,7 @@ static void duplicate_points_mesh(GeometrySet &geometry_set,
   Array<int> offsets = accumulate_counts_to_offsets(selection, counts);
 
   Mesh *new_mesh = BKE_mesh_new_nomain(offsets.last(), 0, 0, 0, 0);
-  MutableSpan<MVert> dst_verts(new_mesh->mvert, new_mesh->totvert);
+  MutableSpan<MVert> dst_verts = new_mesh->vertices_for_write();
 
   threaded_slice_fill(offsets.as_span(), selection, src_verts, dst_verts);
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_curves.cc b/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_curves.cc
index 30b5b7fbd22..ba09acf0bf0 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_curves.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_curves.cc
@@ -3,7 +3,6 @@
 #include "BKE_curves.hh"
 
 #include "DNA_mesh_types.h"
-#include "DNA_meshdata_types.h"
 
 #include "GEO_mesh_to_curve.hh"
 
@@ -23,7 +22,6 @@ static Curves *edge_paths_to_curves_convert(const Mesh &mesh,
                                             const IndexMask start_verts_mask,
                                             const Span<int> next_indices)
 {
-  const Span<MVert> mvert{mesh.mvert, mesh.totvert};
   Vector<int> vert_indices;
   Vector<int> curve_offsets;
   Array<bool> visited(mesh.totvert, false);
diff --git a/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_selection.cc b/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_selection.cc
index 5e9826837a0..ac66e3906a7 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_selection.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_edge_paths_to_selection.cc
@@ -7,9 +7,6 @@
 #include "BLI_set.hh"
 #include "BLI_task.hh"
 
-#include "DNA_mesh_types.h"
-#include "DNA_meshdata_types.h"
-
 #include "node_geometry_util.hh"
 
 #include <set>
@@ -28,6 +25,8 @@ static void edge_paths_to_selection(const Mesh &src_mesh,
                                     const Span<int> next_indices,
                                     MutableSpan<bool> r_selection)
 {
+  const Span<MEdge> edges = src_mesh.edges();
+
   Array<bool> selection(src_mesh.totvert, false);
 
   for (const int start_vert : start_selection) {
@@ -45,8 +44,8 @@ static void edge_paths_to_selection(const Mesh &src_mesh,
     }
   }
 
-  for (const int i : IndexRange(src_mesh.totedge)) {
-    const MEdge &edge = src_mesh.medge[i];
+  for (const int i : edges.index_range()) {
+    const MEdge &edge = edges[i];
     if ((selection[edge.v1] && selection[edge.v2]) &&
         (edge.v1 == next_indices[edge.v2] || edge.v2 == next_indices[edge.v1])) {
       r_selection[i] = true;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_extrude_mesh.cc b/source/blender/nodes/geometry/nodes/node_geo_extrude_mesh.cc
index 237e8ffaa7c..d81b335ebb3 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_extrude_mesh.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_extrude_mesh.cc
@@ -83,31 +83,6 @@ static void save_selection_as_attribute(Mesh &mesh,
   attribute.finish();
 }
 
-static MutableSpan<MVert> mesh_verts(Mesh &mesh)
-{
-  return {mesh.mvert, mesh.totvert};
-}
-static MutableSpan<MEdge> mesh_edges(Mesh &mesh)
-{
-  return {mesh.medge, mesh.totedge};
-}
-static Span<MPoly> mesh_polys(const Mesh &mesh)
-{
-  return {mesh.mpoly, mesh.totpoly};
-}
-static MutableSpan<MPoly> mesh_polys(Mesh &mesh)
-{
-  return {mesh.mpoly, mesh.totpoly};
-}
-static Span<MLoop> mesh_loops(const Mesh &mesh)
-{
-  return {mesh.mloop, mesh.totloop};
-}
-static MutableSpan<MLoop> mesh_loops(Mesh &mesh)
-{
-  return {mesh.mloop, mesh.totloop};
-}
-
 /**
  * \note Some areas in this file rely on the new sections of attributes from #CustomData_realloc
  * to be zeroed.
@@ -142,7 +117,6 @@ static void expand_mesh(Mesh &mesh,
     mesh.totloop += loop_expand;
     CustomData_realloc(&mesh.ldata, mesh.totloop);
   }
-  BKE_mesh_update_customdata_pointers(&mesh, false);
 }
 
 static CustomData &get_customdata(Mesh &mesh, const eAttrDomain domain)
@@ -264,15 +238,15 @@ static void extrude_mesh_vertices(Mesh &mesh,
   const VArray<float3> offsets = evaluator.get_evaluated<float3>(0);
 
   /* This allows parallelizing attribute mixing for new edges. */
-  Array<Vector<int>> vert_to_edge_map = create_vert_to_edge_map(orig_vert_size, mesh_edges(mesh));
+  Array<Vector<int>> vert_to_edge_map = create_vert_to_edge_map(orig_vert_size, mesh.edges());
 
   expand_mesh(mesh, selection.size(), selection.size(), 0, 0);
 
   const IndexRange new_vert_range{orig_vert_size, selection.size()};
   const IndexRange new_edge_range{orig_edge_size, selection.size()};
 
-  MutableSpan<MVert> new_verts = mesh_verts(mesh).slice(new_vert_range);
-  MutableSpan<MEdge> new_edges = mesh_edges(mesh).slice(new_edge_range);
+  MutableSpan<MVert> new_verts = mesh.vertices_for_write().slice(new_vert_range);
+  MutableSpan<MEdge> new_edges = mesh.edges_for_write().slice(new_edge_range);
 
   for (const int i_selection : selection.index_range()) {
     new_edges[i_selection] = new_loose_edge(selection[i_selection], new_vert_range[i_selection]);
@@ -337,8 +311,8 @@ static void extrude_mesh_vertices(Mesh &mesh,
 
 static Array<Vector<int, 2>> mesh_calculate_polys_of_edge(const Mesh &mesh)
 {
-  Span<MPoly> polys = mesh_polys(mesh);
-  Span<MLoop> loops = mesh_loops(mesh);
+  Span<MPoly> polys = mesh.polygons();
+  Span<MLoop> loops = mesh.loops();
   Array<Vector<int, 2>> polys_of_edge(mesh.totedge);
 
   for (const int i_poly : polys.index_range()) {
@@ -396,11 +370,12 @@ template<typename T>
 static VectorSet<int> vert_indices_from_edges(const Mesh &mesh, const Span<T> edge_indices)
 {
   static_assert(is_same_any_v<T, int, int64_t>);
+  const Span<MEdge> edges = mesh.edges();
 
   VectorSet<int> vert_indices;
   vert_indices.reserve(edge_indices.size());
   for (const T i_edge : edge_indices) {
-    const MEdge &edge = mesh.medge[i_edge];
+    const MEdge &edge = edges[i_edge];
     vert_indices.add(edge.v1);
     vert_indices.add(edge.v2);
   }
@@ -413,8 +388,8 @@ static void extrude_mesh_edges(Mesh &mesh,
                                const AttributeOutputs &attribute_outputs)
 {
   const int orig_vert_size = mesh.totvert;
-  Span<MEdge> orig_edges = mesh_edges(mesh);
-  Span<MPoly> orig_polys = mesh_polys(mesh);
+  Span<MEdge> orig_edges = mesh.edges();
+  Span<MPoly> orig_polys = mesh.polygons();
   const int orig_loop_size = mesh.totloop;
 
   bke::MeshFieldContext edge_context{mesh, ATTR_DOMAIN_EDGE};
@@ -463,12 +438,12 @@ static void extrude_mesh_edges(Mesh &mesh,
               new_poly_range.size(),
               new_loop_range.size());
 
-  MutableSpan<MVert> new_verts = mesh_verts(mesh).slice(new_vert_range);
-  MutableSpan<MEdge> connect_edges = mesh_edges(mesh).slice(connect_edge_range);
-  MutableSpan<MEdge> duplicate_edges = mesh_edges(mesh).slice(duplicate_edge_range);
-  MutableSpan<MPoly> polys = mesh_polys(mesh);
+  MutableSpan<MEdge> edges = mesh.edges_for_write();
+  MutableSpan<MEdge> connect_edges = edges.slice(connect_edge_range);
+  MutableSpan<MEdge> duplicate_edges = edges.slice(duplicate_edge_range);
+  MutableSpan<MPoly> polys = mesh.polygons_for_write();
   MutableSpan<MPoly> new_polys = polys.slice(new_poly_range);
-  MutableSpan<MLoop> loops = mesh_loops(mesh);
+  MutableSpan<MLoop> loops = mesh.loops_for_write();
   MutableSpan<MLoop> new_loops = loops.slice(new_loop_range);
 
   for (const int i : connect_edges.index_range()) {
@@ -476,7 +451,7 @@ static void extrude_mesh_edges(Mesh &mesh,
   }
 
   for (const int i : duplicate_edges.index_range()) {
-    const MEdge &orig_edge = mesh.medge[edge_selection[i]];
+    const MEdge &orig_edge = edges[edge_selection[i]];
     const int i_new_vert_1 = new_vert_indices.index_of(orig_edge.v1);
     const int i_new_vert_2 = new_vert_indices.index_of(orig_edge.v2);
     duplicate_edges[i] = new_edge(new_vert_range[i_new_vert_1], new_vert_range[i_new_vert_2]);
@@ -631,6 +606,7 @@ static void extrude_mesh_edges(Mesh &mesh,
     return true;
   });
 
+  MutableSpan<MVert> new_verts = mesh.vertices_for_write().slice(new_vert_range);
   if (edge_offsets.is_single()) {
     const float3 offset = edge_offsets.get_internal_single();
     threading::parallel_for(new_verts.index_range(), 1024, [&](const IndexRange range) {
@@ -676,9 +652,9 @@ static void extrude_mesh_face_regions(Mesh &mesh,
                                       const AttributeOutputs &attribute_outputs)
 {
   const int orig_vert_size = mesh.totvert;
-  Span<MEdge> orig_edges = mesh_edges(mesh);
-  Span<MPoly> orig_polys = mesh_polys(mesh);
-  Span<MLoop> orig_loops = mesh_loops(mesh);
+  Span<MEdge> orig_edges = mesh.edges();
+  Span<MPoly> orig_polys = mesh.polygons();
+  Span<MLoop> orig_loops = mesh.loops();
 
   bke::MeshFieldContext poly_context{mesh, ATTR_DOMAIN_FACE};
   FieldEvaluator poly_evaluator{poly_context, mesh.totpoly};
@@ -781,7 +757,7 @@ static void extrude_mesh_face_regions(Mesh &mesh,
 
   /* The vertices attached to duplicate inner edges also have to be duplicated. */
   for (const int i_edge : new_inner_edge_indices) {
-    const MEdge &edge = mesh.medge[i_edge];
+    const MEdge &edge = orig_edges[i_edge];
     new_vert_indices.add(edge.v1);
     new_vert_indices.add(edge.v2);
   }
@@ -805,13 +781,13 @@ static void extrude_mesh_face_regions(Mesh &mesh,
               side_poly_range.size(),
               side_loop_range.size());
 
-  MutableSpan<MEdge> edges = mesh_edges(mesh);
+  MutableSpan<MEdge> edges = mesh.edges_for_write();
   MutableSpan<MEdge> connect_edges = edges.slice(connect_edge_range);
   MutableSpan<MEdge> boundary_edges = edges.slice(boundary_edge_range);
   MutableSpan<MEdge> new_inner_edges = edges.slice(new_inner_edge_range);
-  MutableSpan<MPoly> polys = mesh_polys(mesh);
+  MutableSpan<MPoly> polys = mesh.polygons_for_write();
   MutableSpan<MPoly> new_polys = polys.slice(side_poly_range);
-  MutableSpan<MLoop> loops = mesh_loops(mesh);
+  MutableSpan<MLoop> loops = mesh.loops_for_write();
   MutableSpan<MLoop> new_loops = loops.slice(side_loop_range);
 
   /* Initialize the edges that form the sides of the extrusion. */
@@ -1000,13 +976,14 @@ static void extrude_mesh_face_regions(Mesh &mesh,
   /* Translate vertices based on the offset. If the vertex is used by a selected edge, it will
    * have been duplicated and only the new vertex should use the offset. Otherwise the vertex might
    * still need an offset, but it was reused on the inside of a region of extruded faces. */
+  MutableSpan<MVert> verts = mesh.vertices_for_write();
   if (poly_offsets.is_single()) {
     const float3 offset = poly_offsets.get_internal_single();
     threading::parallel_for(
         IndexRange(all_selected_verts.size()), 1024, [&](const IndexRange range) {
           for (const int i_orig : all_selected_verts.as_span().slice(range)) {
             const int i_new = new_vert_indices.index_of_try(i_orig);
-            MVert &vert = mesh_verts(mesh)[(i_new == -1) ? i_orig : new_vert_range[i_new]];
+            MVert &vert = verts[(i_new == -1) ? i_orig : new_vert_range[i_new]];
             add_v3_v3(vert.co, offset);
           }
         });
@@ -1017,7 +994,7 @@ static void extrude_mesh_face_regions(Mesh &mesh,
           for (const int i_orig : all_selected_verts.as_span().slice(range)) {
             const int i_new = new_vert_indices.index_of_try(i_orig);
             const float3 offset = vert_offsets[i_orig];
-            MVert &vert = mesh_verts(mesh)[(i_new == -1) ? i_orig : new_vert_range[i_new]];
+            MVert &vert = verts[(i_new == -1) ? i_orig : new_vert_range[i_new]];
             add_v3_v3(vert.co, offset);
           }
         });
@@ -1061,8 +1038,8 @@ static void extrude_individual_mesh_faces(Mesh &mesh,
 {
   const int orig_vert_size = mesh.totvert;
   const int orig_edge_size = mesh.totedge;
-  Span<MPoly> orig_polys = mesh_polys(mesh);
-  Span<MLoop> orig_loops = mesh_loops(mesh);
+  Span<MPoly> orig_polys = mesh.polygons();
+  Span<MLoop> orig_loops = mesh.loops();
 
   /* Use a mesh for the result of the evaluation because the mesh is reallocated before
    * the vertices are moved, and the evaluated result might reference an attribute. */
@@ -1101,13 +1078,13 @@ static void extrude_individual_mesh_faces(Mesh &mesh,
               side_poly_range.size(),
               side_loop_range.size());
 
-  MutableSpan<MVert> new_verts = mesh_verts(mesh).slice(new_vert_range);
-  MutableSpan<MEdge> edges{mesh.medge, mesh.totedge};
+  MutableSpan<MVert> new_verts = mesh.vertices_for_write().slice(new_vert_range);
+  MutableSpan<MEdge> edges = mesh.edges_for_write();
   MutableSpan<MEdge> connect_edges = edges.slice(connect_edge_range);
   MutableSpan<MEdge> duplicate_edges = edges.slice(duplicate_edge_range);
-  MutableSpan<MPoly> polys{mesh.mpoly, mesh.totpoly};
+  MutableSpan<MPoly> polys = mesh.polygons_for_write();
   MutableSpan<MPoly> new_polys = polys.slice(side_poly_range);
-  MutableSpan<MLoop> loops{mesh.mloop, mesh.totloop};
+  MutableSpan<MLoop> loops = mesh.loops_for_write();
 
   /* For every selected polygon, build the faces that form the sides of the extrusion. Filling some
    * of this data like the new edges or polygons could be easily split into separate loops, which
diff --git a/source/blender/nodes/geometry/nodes/node_geo_flip_faces.cc b/source/blender/nodes/geometry/nodes/node_geo_flip_faces.cc
index a752abc2522..2d642ad13d5 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_flip_faces.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_flip_faces.cc
@@ -30,9 +30,8 @@ static void mesh_flip_faces(Mesh &mesh, const Field<bool> &selection_field)
   evaluator.evaluate();
   const IndexMask selection = evaluator.get_evaluated_as_mask(0);
 
-  mesh.mloop = (MLoop *)CustomData_duplicate_referenced_layer(&mesh.ldata, CD_MLOOP, mesh.totloop);
-  const Span<MPoly> polys{mesh.mpoly, mesh.totpoly};
-  MutableSpan<MLoop> loops{mesh.mloop, mesh.totloop};
+  const Span<MPoly> polys = mesh.polygons();
+  MutableSpan<MLoop> loops = mesh.loops_for_write();
 
   for (const int i : selection.index_range()) {
     const MPoly &poly = polys[selection[i]];
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_angle.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_angle.cc
index 3e9fcb10c8e..4b5ea1c8742 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_angle.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_angle.cc
@@ -64,21 +64,20 @@ class AngleFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask UNUSED(mask)) const final
   {
-    Span<MVert> vertices{mesh.mvert, mesh.totvert};
-    Span<MPoly> polys{mesh.mpoly, mesh.totpoly};
-    Span<MLoop> loops{mesh.mloop, mesh.totloop};
+    const Span<MVert> verts = mesh.vertices();
+    const Span<MPoly> polys = mesh.polygons();
+    const Span<MLoop> loops = mesh.loops();
     Array<EdgeMapEntry> edge_map = create_edge_map(polys, loops, mesh.totedge);
 
-    auto angle_fn =
-        [edge_map = std::move(edge_map), vertices, polys, loops](const int i) -> float {
+    auto angle_fn = [edge_map = std::move(edge_map), verts, polys, loops](const int i) -> float {
       if (edge_map[i].face_count != 2) {
         return 0.0f;
       }
       const MPoly &mpoly_1 = polys[edge_map[i].face_index_1];
       const MPoly &mpoly_2 = polys[edge_map[i].face_index_2];
       float3 normal_1, normal_2;
-      BKE_mesh_calc_poly_normal(&mpoly_1, &loops[mpoly_1.loopstart], vertices.data(), normal_1);
-      BKE_mesh_calc_poly_normal(&mpoly_2, &loops[mpoly_2.loopstart], vertices.data(), normal_2);
+      BKE_mesh_calc_poly_normal(&mpoly_1, &loops[mpoly_1.loopstart], verts.data(), normal_1);
+      BKE_mesh_calc_poly_normal(&mpoly_2, &loops[mpoly_2.loopstart], verts.data(), normal_2);
       return angle_normalized_v3v3(normal_1, normal_2);
     };
 
@@ -110,14 +109,14 @@ class SignedAngleFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask UNUSED(mask)) const final
   {
-    const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-    const Span<MEdge> edges(mesh.medge, mesh.totedge);
-    const Span<MPoly> polys(mesh.mpoly, mesh.totpoly);
-    const Span<MLoop> loops(mesh.mloop, mesh.totloop);
+    const Span<MVert> verts = mesh.vertices();
+    const Span<MEdge> edges = mesh.edges();
+    const Span<MPoly> polys = mesh.polygons();
+    const Span<MLoop> loops = mesh.loops();
     Array<EdgeMapEntry> edge_map = create_edge_map(polys, loops, mesh.totedge);
 
     auto angle_fn =
-        [edge_map = std::move(edge_map), vertices, edges, polys, loops](const int i) -> float {
+        [edge_map = std::move(edge_map), verts, edges, polys, loops](const int i) -> float {
       if (edge_map[i].face_count != 2) {
         return 0.0f;
       }
@@ -126,21 +125,18 @@ class SignedAngleFieldInput final : public bke::MeshFieldInput {
 
       /* Find the normals of the 2 polys. */
       float3 poly_1_normal, poly_2_normal;
-      BKE_mesh_calc_poly_normal(
-          &mpoly_1, &loops[mpoly_1.loopstart], vertices.data(), poly_1_normal);
-      BKE_mesh_calc_poly_normal(
-          &mpoly_2, &loops[mpoly_2.loopstart], vertices.data(), poly_2_normal);
+      BKE_mesh_calc_poly_normal(&mpoly_1, &loops[mpoly_1.loopstart], verts.data(), poly_1_normal);
+      BKE_mesh_calc_poly_normal(&mpoly_2, &loops[mpoly_2.loopstart], verts.data(), poly_2_normal);
 
       /* Find the centerpoint of the axis edge */
-      const float3 edge_centerpoint = (float3(vertices[edges[i].v1].co) +
-                                       float3(vertices[edges[i].v2].co)) *
+      const float3 edge_centerpoint = (float3(verts[edges[i].v1].co) +
+                                       float3(verts[edges[i].v2].co)) *
                                       0.5f;
 
       /* Get the centerpoint of poly 2 and subtract the edge centerpoint to get a tangent
        * normal for poly 2. */
       float3 poly_center_2;
-      BKE_mesh_calc_poly_center(
-          &mpoly_2, &loops[mpoly_2.loopstart], vertices.data(), poly_center_2);
+      BKE_mesh_calc_poly_center(&mpoly_2, &loops[mpoly_2.loopstart], verts.data(), poly_center_2);
       const float3 poly_2_tangent = math::normalize(poly_center_2 - edge_centerpoint);
       const float concavity = math::dot(poly_1_normal, poly_2_tangent);
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_neighbors.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_neighbors.cc
index b532b55697b..716cbf589d9 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_neighbors.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_neighbors.cc
@@ -28,9 +28,10 @@ class EdgeNeighborCountFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask UNUSED(mask)) const final
   {
+    const Span<MLoop> loops = mesh.loops();
     Array<int> face_count(mesh.totedge, 0);
-    for (const int i : IndexRange(mesh.totloop)) {
-      face_count[mesh.mloop[i].e]++;
+    for (const MLoop &loop : loops) {
+      face_count[loop.e]++;
     }
 
     return bke::mesh_attributes(mesh).adapt_domain<int>(
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_vertices.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_vertices.cc
index 426e7636d53..7fd2df7c552 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_vertices.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_edge_vertices.cc
@@ -31,7 +31,7 @@ static VArray<int> construct_edge_vertices_gvarray(const Mesh &mesh,
                                                    const VertexNumber vertex,
                                                    const eAttrDomain domain)
 {
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
+  const Span<MEdge> edges = mesh.edges();
   if (domain == ATTR_DOMAIN_EDGE) {
     if (vertex == VERTEX_ONE) {
       return VArray<int>::ForFunc(edges.size(),
@@ -79,19 +79,19 @@ static VArray<float3> construct_edge_positions_gvarray(const Mesh &mesh,
                                                        const VertexNumber vertex,
                                                        const eAttrDomain domain)
 {
-  const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MEdge> edges = mesh.edges();
 
   if (vertex == VERTEX_ONE) {
     return bke::mesh_attributes(mesh).adapt_domain<float3>(
-        VArray<float3>::ForFunc(
-            edges.size(), [vertices, edges](const int i) { return vertices[edges[i].v1].co; }),
+        VArray<float3>::ForFunc(edges.size(),
+                                [verts, edges](const int i) { return verts[edges[i].v1].co; }),
         ATTR_DOMAIN_EDGE,
         domain);
   }
   return bke::mesh_attributes(mesh).adapt_domain<float3>(
       VArray<float3>::ForFunc(edges.size(),
-                              [vertices, edges](const int i) { return vertices[edges[i].v2].co; }),
+                              [verts, edges](const int i) { return verts[edges[i].v2].co; }),
       ATTR_DOMAIN_EDGE,
       domain);
 }
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_area.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_area.cc
index 67b4be0d95d..7dfeaa8e8d9 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_area.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_area.cc
@@ -18,17 +18,17 @@ static void node_declare(NodeDeclarationBuilder &b)
 
 static VArray<float> construct_face_area_varray(const Mesh &mesh, const eAttrDomain domain)
 {
-  const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-  const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
-  const Span<MLoop> loops(mesh.mloop, mesh.totloop);
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
 
-  auto area_fn = [vertices, polygons, loops](const int i) -> float {
-    const MPoly &poly = polygons[i];
-    return BKE_mesh_calc_poly_area(&poly, &loops[poly.loopstart], vertices.data());
+  auto area_fn = [verts, polys, loops](const int i) -> float {
+    const MPoly &poly = polys[i];
+    return BKE_mesh_calc_poly_area(&poly, &loops[poly.loopstart], verts.data());
   };
 
   return bke::mesh_attributes(mesh).adapt_domain<float>(
-      VArray<float>::ForFunc(polygons.size(), area_fn), ATTR_DOMAIN_FACE, domain);
+      VArray<float>::ForFunc(polys.size(), area_fn), ATTR_DOMAIN_FACE, domain);
 }
 
 class FaceAreaFieldInput final : public bke::MeshFieldInput {
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_is_planar.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_is_planar.cc
index 57ab1223d44..b316639fd0a 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_is_planar.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_is_planar.cc
@@ -37,31 +37,30 @@ class PlanarFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask /*mask*/) const final
   {
-    const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-    const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
-    const Span<MLoop> loops(mesh.mloop, mesh.totloop);
+    const Span<MVert> verts = mesh.vertices();
+    const Span<MPoly> polys = mesh.polygons();
+    const Span<MLoop> loops = mesh.loops();
+    const Span<float3> poly_normals{(float3 *)BKE_mesh_poly_normals_ensure(&mesh), mesh.totpoly};
 
     bke::MeshFieldContext context{mesh, ATTR_DOMAIN_FACE};
-    fn::FieldEvaluator evaluator{context, polygons.size()};
+    fn::FieldEvaluator evaluator{context, polys.size()};
     evaluator.add(threshold_);
     evaluator.evaluate();
     const VArray<float> thresholds = evaluator.get_evaluated<float>(0);
 
-    Span<float3> poly_normals{(float3 *)BKE_mesh_poly_normals_ensure(&mesh), polygons.size()};
-
-    auto planar_fn = [vertices, polygons, loops, thresholds, poly_normals](const int i) -> bool {
-      const MPoly &poly = polygons[i];
+    auto planar_fn = [verts, polys, loops, thresholds, poly_normals](const int i) -> bool {
+      const MPoly &poly = polys[i];
       if (poly.totloop <= 3) {
         return true;
       }
       const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
-      float3 reference_normal = poly_normals[i];
+      const float3 &reference_normal = poly_normals[i];
 
       float min = FLT_MAX;
       float max = -FLT_MAX;
 
       for (const int i_loop : poly_loops.index_range()) {
-        const float3 vert = vertices[poly_loops[i_loop].v].co;
+        const float3 vert = verts[poly_loops[i_loop].v].co;
         float dot = math::dot(reference_normal, vert);
         if (dot > max) {
           max = dot;
@@ -74,7 +73,7 @@ class PlanarFieldInput final : public bke::MeshFieldInput {
     };
 
     return bke::mesh_attributes(mesh).adapt_domain<bool>(
-        VArray<bool>::ForFunc(polygons.size(), planar_fn), ATTR_DOMAIN_FACE, domain);
+        VArray<bool>::ForFunc(polys.size(), planar_fn), ATTR_DOMAIN_FACE, domain);
   }
 
   uint64_t hash() const override
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_neighbors.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_neighbors.cc
index c4cb81c5fe5..59d30b997a6 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_neighbors.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_face_neighbors.cc
@@ -21,20 +21,19 @@ static void node_declare(NodeDeclarationBuilder &b)
 
 static VArray<int> construct_neighbor_count_varray(const Mesh &mesh, const eAttrDomain domain)
 {
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
-  const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
-  const Span<MLoop> loops(mesh.mloop, mesh.totloop);
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
 
-  Array<int> edge_count(edges.size(), 0);
-  for (const int i : loops.index_range()) {
-    edge_count[loops[i].e]++;
+  Array<int> edge_count(mesh.totedge, 0);
+  for (const MLoop &loop : loops) {
+    edge_count[loop.e]++;
   }
 
-  Array<int> poly_count(polygons.size(), 0);
-  for (const int poly_i : polygons.index_range()) {
-    const MPoly &poly = polygons[poly_i];
+  Array<int> poly_count(polys.size(), 0);
+  for (const int poly_index : polys.index_range()) {
+    const MPoly &poly = polys[poly_index];
     for (const MLoop &loop : loops.slice(poly.loopstart, poly.totloop)) {
-      poly_count[poly_i] += edge_count[loop.e] - 1;
+      poly_count[poly_index] += edge_count[loop.e] - 1;
     }
   }
 
@@ -71,10 +70,10 @@ class FaceNeighborCountFieldInput final : public bke::MeshFieldInput {
 
 static VArray<int> construct_vertex_count_varray(const Mesh &mesh, const eAttrDomain domain)
 {
-  const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
+  const Span<MPoly> polys = mesh.polygons();
   return bke::mesh_attributes(mesh).adapt_domain<int>(
-      VArray<int>::ForFunc(polygons.size(),
-                           [polygons](const int i) -> float { return polygons[i].totloop; }),
+      VArray<int>::ForFunc(polys.size(),
+                           [polys](const int i) -> float { return polys[i].totloop; }),
       ATTR_DOMAIN_FACE,
       domain);
 }
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_island.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_island.cc
index 5752535d149..53cb3d0a19f 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_island.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_island.cc
@@ -33,7 +33,7 @@ class IslandFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask UNUSED(mask)) const final
   {
-    const Span<MEdge> edges(mesh.medge, mesh.totedge);
+    const Span<MEdge> edges = mesh.edges();
 
     DisjointSet islands(mesh.totvert);
     for (const int i : edges.index_range()) {
@@ -74,7 +74,7 @@ class IslandCountFieldInput final : public bke::MeshFieldInput {
                                  const eAttrDomain domain,
                                  IndexMask UNUSED(mask)) const final
   {
-    const Span<MEdge> edges(mesh.medge, mesh.totedge);
+    const Span<MEdge> edges = mesh.edges();
 
     DisjointSet islands(mesh.totvert);
     for (const int i : edges.index_range()) {
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_vertex_neighbors.cc b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_vertex_neighbors.cc
index 244d454b8d1..ab44a6c8515 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_mesh_vertex_neighbors.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_mesh_vertex_neighbors.cc
@@ -22,8 +22,7 @@ static void node_declare(NodeDeclarationBuilder &b)
 
 static VArray<int> construct_vertex_count_gvarray(const Mesh &mesh, const eAttrDomain domain)
 {
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
-
+  const Span<MEdge> edges = mesh.edges();
   if (domain == ATTR_DOMAIN_POINT) {
     Array<int> counts(mesh.totvert, 0);
     for (const int i : edges.index_range()) {
@@ -63,8 +62,7 @@ class VertexCountFieldInput final : public bke::MeshFieldInput {
 
 static VArray<int> construct_face_count_gvarray(const Mesh &mesh, const eAttrDomain domain)
 {
-  const Span<MLoop> loops(mesh.mloop, mesh.totloop);
-
+  const Span<MLoop> loops = mesh.loops();
   if (domain == ATTR_DOMAIN_POINT) {
     Array<int> vertices(mesh.totvert, 0);
     for (const int i : loops.index_range()) {
diff --git a/source/blender/nodes/geometry/nodes/node_geo_input_shortest_edge_paths.cc b/source/blender/nodes/geometry/nodes/node_geo_input_shortest_edge_paths.cc
index 8549bdfa87d..e13edc8f979 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_input_shortest_edge_paths.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_input_shortest_edge_paths.cc
@@ -2,14 +2,12 @@
 
 #include <queue>
 
-#include "BKE_curves.hh"
-
 #include "BLI_map.hh"
 #include "BLI_math_vec_types.hh"
 #include "BLI_set.hh"
+#include "BLI_task.hh"
 
-#include "DNA_mesh_types.h"
-#include "DNA_meshdata_types.h"
+#include "BKE_mesh.h"
 
 #include "node_geometry_util.hh"
 
@@ -30,7 +28,7 @@ struct EdgeVertMap {
 
   EdgeVertMap(const Mesh &mesh)
   {
-    const Span<MEdge> edges{mesh.medge, mesh.totedge};
+    const Span<MEdge> edges = mesh.edges();
     edges_by_vertex_map.reinitialize(mesh.totvert);
     for (const int edge_i : edges.index_range()) {
       const MEdge &edge = edges[edge_i];
@@ -47,8 +45,7 @@ static void shortest_paths(const Mesh &mesh,
                            MutableSpan<int> r_next_index,
                            MutableSpan<float> r_cost)
 {
-  const Span<MVert> verts{mesh.mvert, mesh.totvert};
-  const Span<MEdge> edges{mesh.medge, mesh.totedge};
+  const Span<MEdge> edges = mesh.edges();
   Array<bool> visited(mesh.totvert, false);
 
   std::priority_queue<VertPriority, std::vector<VertPriority>, std::greater<VertPriority>> queue;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_circle.cc b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_circle.cc
index 9e85547315c..ba3871adc76 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_circle.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_circle.cc
@@ -109,10 +109,10 @@ static Mesh *create_circle_mesh(const float radius,
                                    circle_corner_total(fill_type, verts_num),
                                    circle_face_total(fill_type, verts_num));
   BKE_id_material_eval_ensure_default_slot(&mesh->id);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   /* Assign vertex coordinates. */
   const float angle_delta = 2.0f * (M_PI / static_cast<float>(verts_num));
diff --git a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_cone.cc b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_cone.cc
index cb79ef93de9..98ae1ef1275 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_cone.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_cone.cc
@@ -657,7 +657,7 @@ static Mesh *create_vertex_mesh()
 {
   /* Returns a mesh with a single vertex at the origin. */
   Mesh *mesh = BKE_mesh_new_nomain(1, 0, 0, 0, 0);
-  copy_v3_fl3(mesh->mvert[0].co, 0.0f, 0.0f, 0.0f);
+  copy_v3_fl3(mesh->vertices_for_write().first().co, 0.0f, 0.0f, 0.0f);
   return mesh;
 }
 
@@ -689,10 +689,10 @@ Mesh *create_cylinder_or_cone_mesh(const float radius_top,
       config.tot_verts, config.tot_edges, 0, config.tot_corners, config.tot_faces);
   BKE_id_material_eval_ensure_default_slot(&mesh->id);
 
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   calculate_cone_vertices(verts, config);
   calculate_cone_edges(edges, config);
diff --git a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_grid.cc b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_grid.cc
index 9baf0b3171e..656c5988bef 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_grid.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_grid.cc
@@ -49,10 +49,10 @@ Mesh *create_grid_mesh(const int verts_x,
                                    0,
                                    edges_x * edges_y * 4,
                                    edges_x * edges_y);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   {
     const float dx = edges_x == 0 ? 0.0f : size_x / edges_x;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_line.cc b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_line.cc
index a5edc6c4b3f..130fb8ae589 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_line.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_line.cc
@@ -179,15 +179,15 @@ Mesh *create_line_mesh(const float3 start, const float3 delta, const int count)
 
   Mesh *mesh = BKE_mesh_new_nomain(count, count - 1, 0, 0, 0);
   BKE_id_material_eval_ensure_default_slot(&mesh->id);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
+  MutableSpan<MVert> vertices = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
 
   threading::parallel_invoke(
       1024 < count,
       [&]() {
-        threading::parallel_for(verts.index_range(), 4096, [&](IndexRange range) {
+        threading::parallel_for(vertices.index_range(), 4096, [&](IndexRange range) {
           for (const int i : range) {
-            copy_v3_v3(verts[i].co, start + delta * i);
+            copy_v3_v3(vertices[i].co, start + delta * i);
           }
         });
       },
diff --git a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_uv_sphere.cc b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_uv_sphere.cc
index 85facf1e758..48a95bfcb49 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_uv_sphere.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_mesh_primitive_uv_sphere.cc
@@ -301,10 +301,10 @@ static Mesh *create_uv_sphere_mesh(const float radius, const int segments, const
                                    sphere_corner_total(segments, rings),
                                    sphere_face_total(segments, rings));
   BKE_id_material_eval_ensure_default_slot(&mesh->id);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MEdge> edges{mesh->medge, mesh->totedge};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MEdge> edges = mesh->edges_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   threading::parallel_invoke(
       1024 < segments * rings,
diff --git a/source/blender/nodes/geometry/nodes/node_geo_scale_elements.cc b/source/blender/nodes/geometry/nodes/node_geo_scale_elements.cc
index dbcc5d15fd3..fc467a9424d 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_scale_elements.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_scale_elements.cc
@@ -147,14 +147,22 @@ static float4x4 create_single_axis_transform(const float3 &center,
   return transform;
 }
 
-using GetVertexIndicesFn =
-    FunctionRef<void(const Mesh &mesh, int element_index, VectorSet<int> &r_vertex_indices)>;
+using GetVertexIndicesFn = FunctionRef<void(Span<MEdge> edges,
+                                            Span<MPoly> polys,
+                                            Span<MLoop> loops,
+                                            int element_index,
+                                            VectorSet<int> &r_vertex_indices)>;
 
 static void scale_vertex_islands_uniformly(Mesh &mesh,
                                            const Span<ElementIsland> islands,
                                            const UniformScaleParams &params,
                                            const GetVertexIndicesFn get_vertex_indices)
 {
+  MutableSpan<MVert> verts = mesh.vertices_for_write();
+  const Span<MEdge> edges = mesh.edges();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   threading::parallel_for(islands.index_range(), 256, [&](const IndexRange range) {
     for (const int island_index : range) {
       const ElementIsland &island = islands[island_index];
@@ -164,7 +172,7 @@ static void scale_vertex_islands_uniformly(Mesh &mesh,
 
       VectorSet<int> vertex_indices;
       for (const int poly_index : island.element_indices) {
-        get_vertex_indices(mesh, poly_index, vertex_indices);
+        get_vertex_indices(edges, polys, loops, poly_index, vertex_indices);
         center += params.centers[poly_index];
         scale += params.scales[poly_index];
       }
@@ -175,7 +183,7 @@ static void scale_vertex_islands_uniformly(Mesh &mesh,
       center *= f;
 
       for (const int vert_index : vertex_indices) {
-        MVert &vert = mesh.mvert[vert_index];
+        MVert &vert = verts[vert_index];
         const float3 old_position = vert.co;
         const float3 new_position = transform_with_uniform_scale(old_position, center, scale);
         copy_v3_v3(vert.co, new_position);
@@ -191,6 +199,11 @@ static void scale_vertex_islands_on_axis(Mesh &mesh,
                                          const AxisScaleParams &params,
                                          const GetVertexIndicesFn get_vertex_indices)
 {
+  MutableSpan<MVert> verts = mesh.vertices_for_write();
+  const Span<MEdge> edges = mesh.edges();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   threading::parallel_for(islands.index_range(), 256, [&](const IndexRange range) {
     for (const int island_index : range) {
       const ElementIsland &island = islands[island_index];
@@ -201,7 +214,7 @@ static void scale_vertex_islands_on_axis(Mesh &mesh,
 
       VectorSet<int> vertex_indices;
       for (const int poly_index : island.element_indices) {
-        get_vertex_indices(mesh, poly_index, vertex_indices);
+        get_vertex_indices(edges, polys, loops, poly_index, vertex_indices);
         center += params.centers[poly_index];
         scale += params.scales[poly_index];
         axis += params.axis_vectors[poly_index];
@@ -219,7 +232,7 @@ static void scale_vertex_islands_on_axis(Mesh &mesh,
 
       const float4x4 transform = create_single_axis_transform(center, axis, scale);
       for (const int vert_index : vertex_indices) {
-        MVert &vert = mesh.mvert[vert_index];
+        MVert &vert = verts[vert_index];
         const float3 old_position = vert.co;
         const float3 new_position = transform * old_position;
         copy_v3_v3(vert.co, new_position);
@@ -232,11 +245,14 @@ static void scale_vertex_islands_on_axis(Mesh &mesh,
 
 static Vector<ElementIsland> prepare_face_islands(const Mesh &mesh, const IndexMask face_selection)
 {
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   /* Use the disjoint set data structure to determine which vertices have to be scaled together. */
   DisjointSet disjoint_set(mesh.totvert);
   for (const int poly_index : face_selection) {
-    const MPoly &poly = mesh.mpoly[poly_index];
-    const Span<MLoop> poly_loops{mesh.mloop + poly.loopstart, poly.totloop};
+    const MPoly &poly = polys[poly_index];
+    const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
     for (const int loop_index : IndexRange(poly.totloop - 1)) {
       const int v1 = poly_loops[loop_index].v;
       const int v2 = poly_loops[loop_index + 1].v;
@@ -252,8 +268,8 @@ static Vector<ElementIsland> prepare_face_islands(const Mesh &mesh, const IndexM
 
   /* Gather all of the face indices in each island into separate vectors. */
   for (const int poly_index : face_selection) {
-    const MPoly &poly = mesh.mpoly[poly_index];
-    const Span<MLoop> poly_loops{mesh.mloop + poly.loopstart, poly.totloop};
+    const MPoly &poly = polys[poly_index];
+    const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
     const int island_id = disjoint_set.find_root(poly_loops[0].v);
     const int island_index = island_ids.index_of_or_add(island_id);
     if (island_index == islands.size()) {
@@ -266,10 +282,14 @@ static Vector<ElementIsland> prepare_face_islands(const Mesh &mesh, const IndexM
   return islands;
 }
 
-static void get_face_vertices(const Mesh &mesh, int face_index, VectorSet<int> &r_vertex_indices)
+static void get_face_vertices(const Span<MEdge> /*edges*/,
+                              const Span<MPoly> polys,
+                              const Span<MLoop> loops,
+                              int face_index,
+                              VectorSet<int> &r_vertex_indices)
 {
-  const MPoly &poly = mesh.mpoly[face_index];
-  const Span<MLoop> poly_loops{mesh.mloop + poly.loopstart, poly.totloop};
+  const MPoly &poly = polys[face_index];
+  const Span<MLoop> poly_loops = loops.slice(poly.loopstart, poly.totloop);
   for (const MLoop &loop : poly_loops) {
     r_vertex_indices.add(loop.v);
   }
@@ -290,9 +310,6 @@ static AxisScaleParams evaluate_axis_scale_fields(FieldEvaluator &evaluator,
 
 static void scale_faces_on_axis(Mesh &mesh, const AxisScaleFields &fields)
 {
-  mesh.mvert = static_cast<MVert *>(
-      CustomData_duplicate_referenced_layer(&mesh.vdata, CD_MVERT, mesh.totvert));
-
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_FACE};
   FieldEvaluator evaluator{field_context, mesh.totpoly};
   AxisScaleParams params = evaluate_axis_scale_fields(evaluator, fields);
@@ -315,9 +332,6 @@ static UniformScaleParams evaluate_uniform_scale_fields(FieldEvaluator &evaluato
 
 static void scale_faces_uniformly(Mesh &mesh, const UniformScaleFields &fields)
 {
-  mesh.mvert = static_cast<MVert *>(
-      CustomData_duplicate_referenced_layer(&mesh.vdata, CD_MVERT, mesh.totvert));
-
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_FACE};
   FieldEvaluator evaluator{field_context, mesh.totpoly};
   UniformScaleParams params = evaluate_uniform_scale_fields(evaluator, fields);
@@ -328,10 +342,12 @@ static void scale_faces_uniformly(Mesh &mesh, const UniformScaleFields &fields)
 
 static Vector<ElementIsland> prepare_edge_islands(const Mesh &mesh, const IndexMask edge_selection)
 {
+  const Span<MEdge> edges = mesh.edges();
+
   /* Use the disjoint set data structure to determine which vertices have to be scaled together. */
   DisjointSet disjoint_set(mesh.totvert);
   for (const int edge_index : edge_selection) {
-    const MEdge &edge = mesh.medge[edge_index];
+    const MEdge &edge = edges[edge_index];
     disjoint_set.join(edge.v1, edge.v2);
   }
 
@@ -342,7 +358,7 @@ static Vector<ElementIsland> prepare_edge_islands(const Mesh &mesh, const IndexM
 
   /* Gather all of the edge indices in each island into separate vectors. */
   for (const int edge_index : edge_selection) {
-    const MEdge &edge = mesh.medge[edge_index];
+    const MEdge &edge = edges[edge_index];
     const int island_id = disjoint_set.find_root(edge.v1);
     const int island_index = island_ids.index_of_or_add(island_id);
     if (island_index == islands.size()) {
@@ -355,18 +371,19 @@ static Vector<ElementIsland> prepare_edge_islands(const Mesh &mesh, const IndexM
   return islands;
 }
 
-static void get_edge_vertices(const Mesh &mesh, int edge_index, VectorSet<int> &r_vertex_indices)
+static void get_edge_vertices(const Span<MEdge> edges,
+                              const Span<MPoly> /*polygons*/,
+                              const Span<MLoop> /*loops*/,
+                              int edge_index,
+                              VectorSet<int> &r_vertex_indices)
 {
-  const MEdge &edge = mesh.medge[edge_index];
+  const MEdge &edge = edges[edge_index];
   r_vertex_indices.add(edge.v1);
   r_vertex_indices.add(edge.v2);
 }
 
 static void scale_edges_uniformly(Mesh &mesh, const UniformScaleFields &fields)
 {
-  mesh.mvert = static_cast<MVert *>(
-      CustomData_duplicate_referenced_layer(&mesh.vdata, CD_MVERT, mesh.totvert));
-
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_EDGE};
   FieldEvaluator evaluator{field_context, mesh.totedge};
   UniformScaleParams params = evaluate_uniform_scale_fields(evaluator, fields);
@@ -377,9 +394,6 @@ static void scale_edges_uniformly(Mesh &mesh, const UniformScaleFields &fields)
 
 static void scale_edges_on_axis(Mesh &mesh, const AxisScaleFields &fields)
 {
-  mesh.mvert = static_cast<MVert *>(
-      CustomData_duplicate_referenced_layer(&mesh.vdata, CD_MVERT, mesh.totvert));
-
   bke::MeshFieldContext field_context{mesh, ATTR_DOMAIN_EDGE};
   FieldEvaluator evaluator{field_context, mesh.totedge};
   AxisScaleParams params = evaluate_axis_scale_fields(evaluator, fields);
diff --git a/source/blender/nodes/geometry/nodes/node_geo_set_material.cc b/source/blender/nodes/geometry/nodes/node_geo_set_material.cc
index c6a2f89220a..3aee25b0693 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_set_material.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_set_material.cc
@@ -12,6 +12,7 @@
 #include "DNA_volume_types.h"
 
 #include "BKE_material.h"
+#include "BKE_mesh.h"
 
 namespace blender::nodes::node_geo_set_material_cc {
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_set_position.cc b/source/blender/nodes/geometry/nodes/node_geo_set_position.cc
index d9c7c9422eb..ce453a8ef32 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_set_position.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_set_position.cc
@@ -8,6 +8,7 @@
 #include "DNA_meshdata_types.h"
 
 #include "BKE_curves.hh"
+#include "BKE_mesh.h"
 
 #include "node_geometry_util.hh"
 
@@ -35,14 +36,14 @@ static void set_computed_position_and_offset(GeometryComponent &component,
   switch (component.type()) {
     case GEO_COMPONENT_TYPE_MESH: {
       Mesh *mesh = static_cast<MeshComponent &>(component).get_for_write();
-      MutableSpan<MVert> mverts{mesh->mvert, mesh->totvert};
+      MutableSpan<MVert> verts = mesh->vertices_for_write();
       if (in_positions.is_same(positions.varray)) {
         devirtualize_varray(in_offsets, [&](const auto in_offsets) {
           threading::parallel_for(
               selection.index_range(), grain_size, [&](const IndexRange range) {
                 for (const int i : selection.slice(range)) {
                   const float3 offset = in_offsets[i];
-                  add_v3_v3(mverts[i].co, offset);
+                  add_v3_v3(verts[i].co, offset);
                 }
               });
         });
@@ -54,7 +55,7 @@ static void set_computed_position_and_offset(GeometryComponent &component,
                   selection.index_range(), grain_size, [&](const IndexRange range) {
                     for (const int i : selection.slice(range)) {
                       const float3 new_position = in_positions[i] + in_offsets[i];
-                      copy_v3_v3(mverts[i].co, new_position);
+                      copy_v3_v3(verts[i].co, new_position);
                     }
                   });
             });
diff --git a/source/blender/nodes/geometry/nodes/node_geo_transfer_attribute.cc b/source/blender/nodes/geometry/nodes/node_geo_transfer_attribute.cc
index 539a1488f53..6be37a02bd5 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_transfer_attribute.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_transfer_attribute.cc
@@ -10,6 +10,7 @@
 
 #include "BKE_attribute_math.hh"
 #include "BKE_bvhutils.h"
+#include "BKE_mesh.h"
 #include "BKE_mesh_runtime.h"
 #include "BKE_mesh_sample.hh"
 
@@ -263,6 +264,10 @@ static void get_closest_mesh_corners(const Mesh &mesh,
                                      const MutableSpan<float> r_distances_sq,
                                      const MutableSpan<float3> r_positions)
 {
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   BLI_assert(mesh.totloop > 0);
   Array<int> poly_indices(positions.size());
   get_closest_mesh_polygons(mesh, positions, mask, poly_indices, {}, {});
@@ -270,16 +275,16 @@ static void get_closest_mesh_corners(const Mesh &mesh,
   for (const int i : mask) {
     const float3 position = positions[i];
     const int poly_index = poly_indices[i];
-    const MPoly &poly = mesh.mpoly[poly_index];
+    const MPoly &poly = polys[poly_index];
 
     /* Find the closest vertex in the polygon. */
     float min_distance_sq = FLT_MAX;
     const MVert *closest_mvert;
     int closest_loop_index = 0;
     for (const int loop_index : IndexRange(poly.loopstart, poly.totloop)) {
-      const MLoop &loop = mesh.mloop[loop_index];
+      const MLoop &loop = loops[loop_index];
       const int vertex_index = loop.v;
-      const MVert &mvert = mesh.mvert[vertex_index];
+      const MVert &mvert = verts[vertex_index];
       const float distance_sq = math::distance_squared(position, float3(mvert.co));
       if (distance_sq < min_distance_sq) {
         min_distance_sq = distance_sq;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc b/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc
index a59704291cd..5900e234220 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_uv_pack_islands.cc
@@ -5,6 +5,8 @@
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
 
+#include "BKE_mesh.h"
+
 #include "node_geometry_util.hh"
 
 namespace blender::nodes::node_geo_uv_pack_islands_cc {
@@ -35,8 +37,12 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
                                            const float margin,
                                            const eAttrDomain domain)
 {
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   bke::MeshFieldContext face_context{mesh, ATTR_DOMAIN_FACE};
-  FieldEvaluator face_evaluator{face_context, mesh.totpoly};
+  FieldEvaluator face_evaluator{face_context, polys.size()};
   face_evaluator.add(selection_field);
   face_evaluator.evaluate();
   const IndexMask selection = face_evaluator.get_evaluated_as_mask(0);
@@ -50,14 +56,9 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
   evaluator.add_with_destination(uv_field, uv.as_mutable_span());
   evaluator.evaluate();
 
-  const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
-  const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
-  const Span<MLoop> loops(mesh.mloop, mesh.totloop);
-
   ParamHandle *handle = GEO_uv_parametrizer_construct_begin();
   for (const int mp_index : selection) {
-    const MPoly &mp = polygons[mp_index];
+    const MPoly &mp = polys[mp_index];
     Array<ParamKey, 16> mp_vkeys(mp.totloop);
     Array<bool, 16> mp_pin(mp.totloop);
     Array<bool, 16> mp_select(mp.totloop);
@@ -66,7 +67,7 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
     for (const int i : IndexRange(mp.totloop)) {
       const MLoop &ml = loops[mp.loopstart + i];
       mp_vkeys[i] = ml.v;
-      mp_co[i] = vertices[ml.v].co;
+      mp_co[i] = verts[ml.v].co;
       mp_uv[i] = uv[mp.loopstart + i];
       mp_pin[i] = false;
       mp_select[i] = false;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc b/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc
index 786438ad62a..3095cac6a50 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_uv_unwrap.cc
@@ -5,6 +5,8 @@
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
 
+#include "BKE_mesh.h"
+
 #include "UI_interface.h"
 #include "UI_resources.h"
 
@@ -60,8 +62,13 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
                                            const GeometryNodeUVUnwrapMethod method,
                                            const eAttrDomain domain)
 {
+  const Span<MVert> verts = mesh.vertices();
+  const Span<MEdge> edges = mesh.edges();
+  const Span<MPoly> polys = mesh.polygons();
+  const Span<MLoop> loops = mesh.loops();
+
   bke::MeshFieldContext face_context{mesh, ATTR_DOMAIN_FACE};
-  FieldEvaluator face_evaluator{face_context, mesh.totpoly};
+  FieldEvaluator face_evaluator{face_context, polys.size()};
   face_evaluator.add(selection_field);
   face_evaluator.evaluate();
   const IndexMask selection = face_evaluator.get_evaluated_as_mask(0);
@@ -70,21 +77,16 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
   }
 
   bke::MeshFieldContext edge_context{mesh, ATTR_DOMAIN_EDGE};
-  FieldEvaluator edge_evaluator{edge_context, mesh.totedge};
+  FieldEvaluator edge_evaluator{edge_context, edges.size()};
   edge_evaluator.add(seam_field);
   edge_evaluator.evaluate();
   const IndexMask seam = edge_evaluator.get_evaluated_as_mask(0);
 
-  const Span<MVert> vertices(mesh.mvert, mesh.totvert);
-  const Span<MEdge> edges(mesh.medge, mesh.totedge);
-  const Span<MPoly> polygons(mesh.mpoly, mesh.totpoly);
-  const Span<MLoop> loops(mesh.mloop, mesh.totloop);
-
   Array<float3> uv(loops.size(), float3(0));
 
   ParamHandle *handle = GEO_uv_parametrizer_construct_begin();
   for (const int mp_index : selection) {
-    const MPoly &mp = polygons[mp_index];
+    const MPoly &mp = polys[mp_index];
     Array<ParamKey, 16> mp_vkeys(mp.totloop);
     Array<bool, 16> mp_pin(mp.totloop);
     Array<bool, 16> mp_select(mp.totloop);
@@ -93,7 +95,7 @@ static VArray<float3> construct_uv_gvarray(const Mesh &mesh,
     for (const int i : IndexRange(mp.totloop)) {
       const MLoop &ml = loops[mp.loopstart + i];
       mp_vkeys[i] = ml.v;
-      mp_co[i] = vertices[ml.v].co;
+      mp_co[i] = verts[ml.v].co;
       mp_uv[i] = uv[mp.loopstart + i];
       mp_pin[i] = false;
       mp_select[i] = false;
diff --git a/source/blender/nodes/geometry/nodes/node_geo_volume_to_mesh.cc b/source/blender/nodes/geometry/nodes/node_geo_volume_to_mesh.cc
index 91429560ac8..5788a744041 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_volume_to_mesh.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_volume_to_mesh.cc
@@ -123,9 +123,9 @@ static Mesh *create_mesh_from_volume_grids(Span<openvdb::GridBase::ConstPtr> gri
 
   Mesh *mesh = BKE_mesh_new_nomain(vert_offset, 0, 0, loop_offset, poly_offset);
   BKE_id_material_eval_ensure_default_slot(&mesh->id);
-  MutableSpan<MVert> verts{mesh->mvert, mesh->totvert};
-  MutableSpan<MLoop> loops{mesh->mloop, mesh->totloop};
-  MutableSpan<MPoly> polys{mesh->mpoly, mesh->totpoly};
+  MutableSpan<MVert> verts = mesh->vertices_for_write();
+  MutableSpan<MPoly> polys = mesh->polygons_for_write();
+  MutableSpan<MLoop> loops = mesh->loops_for_write();
 
   for (const int i : grids.index_range()) {
     const bke::OpenVDBMeshData &data = mesh_data[i];