From 7499f96805dc3c85e3d33d08c80146490267ceaa Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Mon, 19 Dec 2022 12:41:47 -0500
Subject: [PATCH] Improve memory usage when reading ExodusII grids.

For my big heart mesh (which has about 2.4 million elements) this lowers
peak memory consumption from 1.3 GB to 1.2 GB.
---
 include/deal.II/grid/connectivity.h |  2 +-
 source/grid/grid_in.cc              | 59 +++++++++++++++++------------
 2 files changed, 35 insertions(+), 26 deletions(-)

diff --git a/include/deal.II/grid/connectivity.h b/include/deal.II/grid/connectivity.h
index d4264bf045..669a5b41fe 100644
--- a/include/deal.II/grid/connectivity.h
+++ b/include/deal.II/grid/connectivity.h
@@ -1130,7 +1130,7 @@ namespace internal
 
           std::sort(keys.begin(), keys.end());
 
-          ptr_0.reserve(n_unique_entities);
+          ptr_0.reserve(n_unique_entities + 1);
           col_0.reserve(n_unique_entity_vertices);
         }
 
diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index cd75e69224..312400bc28 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -3616,6 +3616,7 @@ namespace
         // Collect into a sortable data structure:
         std::vector<std::pair<std::size_t, std::vector<int>>>
           face_id_to_side_sets;
+        face_id_to_side_sets.reserve(face_side_sets.size());
         for (auto &pair : face_side_sets)
           {
             Assert(pair.second.size() > 0, ExcInternalError());
@@ -3632,6 +3633,10 @@ namespace
                                                         b.second.end());
                   });
 
+        if (dim == 2)
+          subcelldata.boundary_lines.reserve(face_id_to_side_sets.size());
+        else if (dim == 3)
+          subcelldata.boundary_quads.reserve(face_id_to_side_sets.size());
         types::boundary_id current_b_or_m_id = 0;
         for (const auto &pair : face_id_to_side_sets)
           {
@@ -3741,13 +3746,7 @@ GridIn<dim, spacedim>::read_exodusii(
   AssertDimension(mesh_dimension, spacedim);
 
   // Read nodes:
-  std::vector<double> xs(n_nodes);
-  std::vector<double> ys(n_nodes);
-  std::vector<double> zs(n_nodes);
-
-  ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
-  AssertThrowExodusII(ierr);
-
+  //
   // Even if there is a node numbering array the values stored inside the
   // ExodusII file must use the contiguous, internal ordering (see Section 4.5
   // of the manual - "Internal (contiguously numbered) node and element IDs
@@ -3756,29 +3755,39 @@ GridIn<dim, spacedim>::read_exodusii(
   // connectivity.")
   std::vector<Point<spacedim>> vertices;
   vertices.reserve(n_nodes);
-  for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
-    {
-      switch (spacedim)
-        {
-          case 1:
-            vertices.emplace_back(xs[vertex_n]);
-            break;
-          case 2:
-            vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
-            break;
-          case 3:
-            vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
-            break;
-          default:
-            Assert(spacedim <= 3, ExcNotImplemented());
-        }
-    }
+  {
+    std::vector<double> xs(n_nodes);
+    std::vector<double> ys(n_nodes);
+    std::vector<double> zs(n_nodes);
+
+    ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
+    AssertThrowExodusII(ierr);
+
+    for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
+      {
+        switch (spacedim)
+          {
+            case 1:
+              vertices.emplace_back(xs[vertex_n]);
+              break;
+            case 2:
+              vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
+              break;
+            case 3:
+              vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
+              break;
+            default:
+              Assert(spacedim <= 3, ExcNotImplemented());
+          }
+      }
+  }
 
   std::vector<int> element_block_ids(n_element_blocks);
   ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
   AssertThrowExodusII(ierr);
 
   std::vector<CellData<dim>> cells;
+  cells.reserve(n_elements);
   // Elements are grouped together by same reference cell type in element
   // blocks. There may be multiple blocks for a single reference cell type,
   // but "each element block may contain only one element type".
@@ -3829,7 +3838,7 @@ GridIn<dim, spacedim>::read_exodusii(
                 connection[elem_n + i] - 1;
             }
           cell.material_id = element_block_id;
-          cells.push_back(cell);
+          cells.push_back(std::move(cell));
         }
     }
 
-- 
2.39.5