diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C
index 4667a64470f..c83aa41db71 100644
--- a/src/mesh/mesh_tools.C
+++ b/src/mesh/mesh_tools.C
@@ -70,8 +70,8 @@ public:
 
   void operator()(const ConstElemRange & range)
   {
-    for (ConstElemRange::const_iterator it = range.begin(); it !=range.end(); ++it)
-      _weight += (*it)->n_nodes();
+    for (const auto & elem : range)
+      _weight += elem->n_nodes();
   }
 
   dof_id_type weight() const
@@ -107,22 +107,18 @@ public:
 
   void operator()(const ConstNodeRange & range)
   {
-    for (ConstNodeRange::const_iterator it = range.begin(); it != range.end(); ++it)
+    for (const auto & node : range)
       {
-        const Node * node = *it;
         libmesh_assert(node);
-
         _bbox.union_with(*node);
       }
   }
 
   void operator()(const ConstElemRange & range)
   {
-    for (ConstElemRange::const_iterator it = range.begin(); it != range.end(); ++it)
+    for (const auto & elem : range)
       {
-        const Elem * elem = *it;
         libmesh_assert(elem);
-
         _bbox.union_with(elem->loose_bounding_box());
       }
   }
@@ -540,12 +536,11 @@ dof_id_type MeshTools::n_non_subactive_elem_of_type_at_level(const MeshBase & me
                                                              const unsigned int level)
 {
   dof_id_type cnt = 0;
-  // iterate over the elements of the specified type
-  MeshBase::const_element_iterator el = mesh.type_elements_begin(type);
-  const MeshBase::const_element_iterator end = mesh.type_elements_end(type);
 
-  for (; el!=end; ++el)
-    if (((*el)->level() == level) && !(*el)->subactive())
+  // iterate over the elements of the specified type
+  for (const auto & elem : as_range(std::make_pair(mesh.type_elements_begin(type),
+                                                   mesh.type_elements_end(type))))
+    if ((elem->level() == level) && !elem->subactive())
       cnt++;
 
   return cnt;
diff --git a/src/mesh/nemesis_io_helper.C b/src/mesh/nemesis_io_helper.C
index 5f521100a6c..e5eb80286ac 100644
--- a/src/mesh/nemesis_io_helper.C
+++ b/src/mesh/nemesis_io_helper.C
@@ -1350,12 +1350,8 @@ Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(const Me
   if (verbose)
     {
       libMesh::out << "[" << this->processor_id() << "] neighboring_processor_ids = ";
-      for (std::set<unsigned>::iterator it = neighboring_processor_ids.begin();
-           it != neighboring_processor_ids.end();
-           ++it)
-        {
-          libMesh::out << *it << " ";
-        }
+      for (const auto & id : neighboring_processor_ids)
+        libMesh::out << id << " ";
       libMesh::out << std::endl;
     }
 
@@ -1372,13 +1368,12 @@ Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(const Me
   if (verbose)
     {
       // Print out counts of border elements for each processor
-      for (proc_border_elem_sets_iterator it=this->proc_border_elem_sets.begin();
-           it != this->proc_border_elem_sets.end(); ++it)
+      for (const auto & pr : proc_border_elem_sets)
         {
           libMesh::out << "[" << this->processor_id() << "] "
                        << "Proc "
-                       << it->first << " communicates "
-                       << it->second.size() << " elements." << std::endl;
+                       << pr.first << " communicates "
+                       << pr.second.size() << " elements." << std::endl;
         }
     }
 
@@ -1544,10 +1539,8 @@ void Nemesis_IO_Helper::compute_num_global_nodesets(const MeshBase & pmesh)
       libMesh::out << std::endl;
 
       libMesh::out << "[" << this->processor_id() << "] local_node_boundary_ids = ";
-      for (std::set<boundary_id_type>::iterator it = local_node_boundary_ids.begin();
-           it != local_node_boundary_ids.end();
-           ++it)
-        libMesh::out << *it << ", ";
+      for (const auto & id : local_node_boundary_ids)
+        libMesh::out << id << ", ";
       libMesh::out << std::endl;
     }
 
@@ -1683,13 +1676,11 @@ void Nemesis_IO_Helper::compute_num_global_elem_blocks(const MeshBase & pmesh)
   if (verbose)
     {
       libMesh::out << "[" << this->processor_id() << "] ";
-      for (std::map<subdomain_id_type, unsigned>::iterator it=global_subdomain_counts.begin();
-           it != global_subdomain_counts.end();
-           ++it)
+      for (const auto & pr : global_subdomain_counts)
         {
           libMesh::out << "ID: "
-                       << static_cast<unsigned>(it->first)
-                       << ", Count: " << it->second << ", ";
+                       << static_cast<unsigned>(pr.first)
+                       << ", Count: " << pr.second << ", ";
         }
       libMesh::out << std::endl;
     }
@@ -1700,12 +1691,9 @@ void Nemesis_IO_Helper::compute_num_global_elem_blocks(const MeshBase & pmesh)
   this->global_elem_blk_cnts.resize(global_subdomain_ids.size());
 
   unsigned cnt=0;
-  for (std::set<subdomain_id_type>::iterator it=global_subdomain_ids.begin();
-       it != global_subdomain_ids.end(); ++it)
-    {
-      // Find the entry in the local map, note: if not found, will be created with 0 default value, which is OK...
-      this->global_elem_blk_cnts[cnt++] = global_subdomain_counts[*it];
-    }
+  // Find the entry in the local map, note: if not found, will be created with 0 default value, which is OK...
+  for (const auto & id : global_subdomain_ids)
+    this->global_elem_blk_cnts[cnt++] = global_subdomain_counts[id];
 
   // Sum up subdomain counts from all processors
   this->comm().sum(this->global_elem_blk_cnts);
@@ -1760,29 +1748,18 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh)
   // Make sure there is no leftover information in the subdomain_map, and reserve
   // enough space to store the elements we need.
   this->subdomain_map.clear();
-  for (std::map<subdomain_id_type, unsigned>::iterator it=this->local_subdomain_counts.begin();
-       it != this->local_subdomain_counts.end();
-       ++it)
+  for (const auto & pr : local_subdomain_counts)
     {
-      subdomain_id_type cur_subdomain = it->first;
-
-      /*
-      // We can't have a zero subdomain ID in Exodus (for some reason?)
-      // so map zero subdomains to a max value...
-      if (cur_subdomain == 0)
-      cur_subdomain = std::numeric_limits<subdomain_id_type>::max();
-      */
-
       if (verbose)
         {
           libMesh::out << "[" << this->processor_id() << "] "
-                       << "local_subdomain_counts [" << static_cast<unsigned>(cur_subdomain) << "]= "
-                       << it->second
+                       << "local_subdomain_counts [" << static_cast<unsigned>(pr.first) << "]= "
+                       << pr.second
                        << std::endl;
         }
 
-      // *it.first is the subdomain ID, *it.second is the number of elements it contains
-      this->subdomain_map[ cur_subdomain ].reserve( it->second );
+      // pr.first is the subdomain ID, pr.second is the number of elements it contains
+      this->subdomain_map[pr.first].reserve(pr.second);
     }
 
 
@@ -1812,20 +1789,18 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh)
   this->libmesh_node_num_to_exodus.clear();
 
   // Set the map for nodes
-  for (std::set<int>::iterator it = this->nodes_attached_to_local_elems.begin();
-       it != this->nodes_attached_to_local_elems.end();
-       ++it)
+  for (const auto & id : nodes_attached_to_local_elems)
     {
       // I.e. given exodus_node_id,
       // exodus_node_num_to_libmesh[ exodus_node_id ] returns the libmesh ID for that node.
       // Note that even though most of Exodus is 1-based, this code will map an Exodus ID of
       // zero to some libmesh node ID.  Is that a problem?
-      this->exodus_node_num_to_libmesh.push_back(*it);
+      this->exodus_node_num_to_libmesh.push_back(id);
 
       // Likewise, given libmesh_node_id,
       // libmesh_node_num_to_exodus[ libmesh_node_id ] returns the *Exodus* ID for that node.
       // Unlike the exodus_node_num_to_libmesh vector above, this one is a std::map
-      this->libmesh_node_num_to_exodus[*it] =
+      this->libmesh_node_num_to_exodus[id] =
         cast_int<int>(this->exodus_node_num_to_libmesh.size()); // should never be zero...
     }
 
@@ -1838,18 +1813,16 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh)
   this->libmesh_elem_num_to_exodus.clear();
 
   // Now loop over each subdomain and get a unique numbering for the elements
-  for (std::map<subdomain_id_type, std::vector<unsigned int>>::iterator it = this->subdomain_map.begin();
-       it != this->subdomain_map.end();
-       ++it)
+  for (auto & pr : subdomain_map)
     {
-      block_ids.push_back(it->first);
+      block_ids.push_back(pr.first);
 
       // Vector of element IDs for this subdomain
-      std::vector<unsigned int> & elem_ids_this_subdomain = it->second;
+      std::vector<unsigned int> & elem_ids_this_subdomain = pr.second;
 
       // The code below assumes this subdomain block is not empty, make sure that's the case!
       if (elem_ids_this_subdomain.size() == 0)
-        libmesh_error_msg("Error, no element IDs found in subdomain " << it->first);
+        libmesh_error_msg("Error, no element IDs found in subdomain " << pr.first);
 
       ExodusII_IO_Helper::ElementMaps em;
 
@@ -1863,7 +1836,7 @@ void Nemesis_IO_Helper::build_element_and_node_maps(const MeshBase & pmesh)
 
       // Get a reference to the connectivity vector for this subdomain.  This vector
       // is most likely empty, we are going to fill it up now.
-      std::vector<int> & current_block_connectivity = this->block_id_to_elem_connectivity[it->first];
+      std::vector<int> & current_block_connectivity = this->block_id_to_elem_connectivity[pr.first];
 
       // Just in case it's not already empty...
       current_block_connectivity.clear();
@@ -1961,33 +1934,27 @@ void Nemesis_IO_Helper::compute_border_node_ids(const MeshBase & pmesh)
                      << " sets of nodes."
                      << std::endl;
 
-        for (proc_nodes_touched_iterator it = proc_nodes_touched.begin();
-             it != proc_nodes_touched.end();
-             ++it)
-          {
-            libMesh::out << "[" << this->processor_id()
-                         << "] proc_nodes_touched[" << it->first << "] has "
-                         << it->second.size()
-                         << " entries."
-                         << std::endl;
-          }
+        for (const auto & pr : proc_nodes_touched)
+          libMesh::out << "[" << this->processor_id()
+                       << "] proc_nodes_touched[" << pr.first << "] has "
+                       << pr.second.size()
+                       << " entries."
+                       << std::endl;
       }
 
 
     // Loop over all the sets we just created and compute intersections with the
     // this processor's set.  Obviously, don't intersect with ourself.
-    for (proc_nodes_touched_iterator it = proc_nodes_touched.begin();
-         it != proc_nodes_touched.end();
-         ++it)
+    for (auto & pr : proc_nodes_touched)
       {
         // Don't compute intersections with ourself
-        if (it->first == this->processor_id())
+        if (pr.first == this->processor_id())
           continue;
 
         // Otherwise, compute intersection with other processor and ourself
         std::set<unsigned> & my_set = proc_nodes_touched[this->processor_id()];
-        std::set<unsigned> & other_set = it->second;
-        std::set<unsigned> & result_set = this->proc_nodes_touched_intersections[ it->first ]; // created if does not exist
+        std::set<unsigned> & other_set = pr.second;
+        std::set<unsigned> & result_set = this->proc_nodes_touched_intersections[pr.first]; // created if does not exist
 
         std::set_intersection(my_set.begin(), my_set.end(),
                               other_set.begin(), other_set.end(),
@@ -1996,25 +1963,19 @@ void Nemesis_IO_Helper::compute_border_node_ids(const MeshBase & pmesh)
 
     if (verbose)
       {
-        for (proc_nodes_touched_iterator it = this->proc_nodes_touched_intersections.begin();
-             it != this->proc_nodes_touched_intersections.end();
-             ++it)
-          {
-            libMesh::out << "[" << this->processor_id()
-                         << "] this->proc_nodes_touched_intersections[" << it->first << "] has "
-                         << it->second.size()
-                         << " entries."
-                         << std::endl;
-          }
+        for (const auto & pr : proc_nodes_touched_intersections)
+          libMesh::out << "[" << this->processor_id()
+                       << "] this->proc_nodes_touched_intersections[" << pr.first << "] has "
+                       << pr.second.size()
+                       << " entries."
+                       << std::endl;
       }
 
     // Compute the set_union of all the preceding intersections.  This will be the set of
     // border node IDs for this processor.
-    for (proc_nodes_touched_iterator it = this->proc_nodes_touched_intersections.begin();
-         it != this->proc_nodes_touched_intersections.end();
-         ++it)
+    for (auto & pr : proc_nodes_touched_intersections)
       {
-        std::set<unsigned> & other_set = it->second;
+        std::set<unsigned> & other_set = pr.second;
         std::set<unsigned> intermediate_result; // Don't think we can insert into one of the sets we're unioning...
 
         std::set_union(this->border_node_ids.begin(), this->border_node_ids.end(),
@@ -2096,19 +2057,13 @@ void Nemesis_IO_Helper::write_nodesets(const MeshBase & mesh)
   // See what we got
   if (verbose)
     {
-      for (std::map<boundary_id_type, std::vector<int>>::iterator it = local_node_boundary_id_lists.begin();
-           it != local_node_boundary_id_lists.end();
-           ++it)
+      for (const auto & pr : local_node_boundary_id_lists)
         {
-          libMesh::out << "[" << this->processor_id() << "] ID: " << it->first << ", ";
-
-          std::vector<int> & current_id_set = it->second;
+          libMesh::out << "[" << this->processor_id() << "] ID: " << pr.first << ", ";
 
           // Libmesh node ID (Exodus Node ID)
-          for (std::size_t j=0; j<current_id_set.size(); ++j)
-            libMesh::out << current_id_set[j]
-                         << ", ";
-
+          for (const auto & id : pr.second)
+            libMesh::out << id << ", ";
           libMesh::out << std::endl;
         }
     }
diff --git a/src/mesh/replicated_mesh.C b/src/mesh/replicated_mesh.C
index c7a01a47e72..cee389dff27 100644
--- a/src/mesh/replicated_mesh.C
+++ b/src/mesh/replicated_mesh.C
@@ -575,34 +575,23 @@ void ReplicatedMesh::clear ()
   // Call parent clear function
   MeshBase::clear();
 
-
   // Clear our elements and nodes
-  {
-    std::vector<Elem *>::iterator       it  = _elements.begin();
-    const std::vector<Elem *>::iterator end = _elements.end();
+  // There is no need to remove the elements from
+  // the BoundaryInfo data structure since we
+  // already cleared it.
+  for (auto & elem : _elements)
+    delete elem;
 
-    // There is no need to remove the elements from
-    // the BoundaryInfo data structure since we
-    // already cleared it.
-    for (; it != end; ++it)
-      delete *it;
-
-    _elements.clear();
-  }
+  _elements.clear();
 
   // clear the nodes data structure
-  {
-    std::vector<Node *>::iterator       it  = _nodes.begin();
-    const std::vector<Node *>::iterator end = _nodes.end();
-
-    // There is no need to remove the nodes from
-    // the BoundaryInfo data structure since we
-    // already cleared it.
-    for (; it != end; ++it)
-      delete *it;
+  // There is no need to remove the nodes from
+  // the BoundaryInfo data structure since we
+  // already cleared it.
+  for (auto & node : _nodes)
+    delete node;
 
-    _nodes.clear();
-  }
+  _nodes.clear();
 }
 
 
@@ -751,7 +740,6 @@ void ReplicatedMesh::renumber_nodes_and_elements ()
       // are not connected to any elements and may be deleted
       // if desired.
 
-      // (This code block will erase the unused nodes)
       // Now, delete the unused nodes
       {
         std::vector<Node *>::iterator nd        = _nodes.begin();
@@ -759,21 +747,20 @@ void ReplicatedMesh::renumber_nodes_and_elements ()
 
         std::advance (nd, next_free_node);
 
-        for (std::vector<Node *>::iterator it=nd;
-             it != end; ++it)
+        for (auto & node : as_range(std::make_pair(nd, end)))
           {
             // Mesh modification code might have already deleted some
             // nodes
-            if (*it == libmesh_nullptr)
+            if (node == libmesh_nullptr)
               continue;
 
             // remove any boundary information associated with
             // this node
-            this->get_boundary_info().remove (*it);
+            this->get_boundary_info().remove (node);
 
             // delete the node
-            delete *it;
-            *it = libmesh_nullptr;
+            delete node;
+            node = libmesh_nullptr;
           }
 
         _nodes.erase (nd, end);
@@ -922,12 +909,8 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
             // Container to catch boundary IDs passed back from BoundaryInfo.
             std::vector<boundary_id_type> bc_ids;
 
-            MeshBase::element_iterator elem_it  = mesh_array[i]->elements_begin();
-            MeshBase::element_iterator elem_end = mesh_array[i]->elements_end();
-            for ( ; elem_it != elem_end; ++elem_it)
+            for (auto & el : mesh_array[i]->element_ptr_range())
               {
-                Elem * el = *elem_it;
-
                 // Now check whether elem has a face on the specified boundary
                 for (auto side_id : el->side_index_range())
                   if (el->neighbor_ptr(side_id) == libmesh_nullptr)
@@ -1084,21 +1067,14 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
           // Otherwise, use a simple N^2 search to find the closest matching points. This can be helpful
           // in the case that we have tolerance issues which cause mismatch between the two surfaces
           // that are being stitched.
-
-          std::set<dof_id_type>::iterator set_it     = this_boundary_node_ids.begin();
-          std::set<dof_id_type>::iterator set_it_end = this_boundary_node_ids.end();
-          for ( ; set_it != set_it_end; ++set_it)
+          for (const auto & this_node_id : this_boundary_node_ids)
             {
-              dof_id_type this_node_id = *set_it;
               Node & this_node = this->node_ref(this_node_id);
 
               bool found_matching_nodes = false;
 
-              std::set<dof_id_type>::iterator other_set_it     = other_boundary_node_ids.begin();
-              std::set<dof_id_type>::iterator other_set_it_end = other_boundary_node_ids.end();
-              for ( ; other_set_it != other_set_it_end; ++other_set_it)
+              for (const auto & other_node_id : other_boundary_node_ids)
                 {
-                  dof_id_type other_node_id = *other_set_it;
                   Node & other_node = other_mesh->node_ref(other_node_id);
 
                   Real node_distance = (this_node - other_node).norm();
@@ -1119,32 +1095,26 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
         }
 
       // Build up the node_to_elems_map, using only one loop over other_mesh
-      {
-        MeshBase::element_iterator other_elem_it  = other_mesh->elements_begin();
-        MeshBase::element_iterator other_elem_end = other_mesh->elements_end();
-        for (; other_elem_it != other_elem_end; ++other_elem_it)
-          {
-            Elem * el = *other_elem_it;
-
-            // For each node on the element, find the corresponding node
-            // on "this" Mesh, 'this_node_id', if it exists, and push
-            // the current element ID back onto node_to_elems_map[this_node_id].
-            // For that we will use the reverse mapping we created at
-            // the same time as the forward mapping.
-            for (auto & n : el->node_ref_range())
-              {
-                dof_id_type other_node_id = n.id();
-                std::map<dof_id_type, dof_id_type>::iterator it =
-                  other_to_this_node_map.find(other_node_id);
+      for (auto & el : other_mesh->element_ptr_range())
+        {
+          // For each node on the element, find the corresponding node
+          // on "this" Mesh, 'this_node_id', if it exists, and push
+          // the current element ID back onto node_to_elems_map[this_node_id].
+          // For that we will use the reverse mapping we created at
+          // the same time as the forward mapping.
+          for (auto & n : el->node_ref_range())
+            {
+              dof_id_type other_node_id = n.id();
+              std::map<dof_id_type, dof_id_type>::iterator it =
+                other_to_this_node_map.find(other_node_id);
 
-                if (it != other_to_this_node_map.end())
-                  {
-                    dof_id_type this_node_id = it->second;
-                    node_to_elems_map[this_node_id].push_back( el->id() );
-                  }
-              }
-          }
-      }
+              if (it != other_to_this_node_map.end())
+                {
+                  dof_id_type this_node_id = it->second;
+                  node_to_elems_map[this_node_id].push_back( el->id() );
+                }
+            }
+        }
 
       if (verbose)
         {
@@ -1183,52 +1153,33 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
       LOG_SCOPE("stitch_meshes copying", "ReplicatedMesh");
 
       // need to increment node and element IDs of other_mesh before copying to this mesh
-      MeshBase::node_iterator node_it  = other_mesh->nodes_begin();
-      MeshBase::node_iterator node_end = other_mesh->nodes_end();
-      for (; node_it != node_end; ++node_it)
+      for (auto & nd : other_mesh->node_ptr_range())
         {
-          Node * nd = *node_it;
           dof_id_type new_id = nd->id() + node_delta;
           nd->set_id(new_id);
         }
 
-      MeshBase::element_iterator elem_it  = other_mesh->elements_begin();
-      MeshBase::element_iterator elem_end = other_mesh->elements_end();
-      for (; elem_it != elem_end; ++elem_it)
+      for (auto & el : other_mesh->element_ptr_range())
         {
-          Elem * el = *elem_it;
           dof_id_type new_id = el->id() + elem_delta;
           el->set_id(new_id);
         }
 
       // Also, increment the node_to_node_map and node_to_elems_map
-      std::map<dof_id_type, dof_id_type>::iterator node_map_it     = node_to_node_map.begin();
-      std::map<dof_id_type, dof_id_type>::iterator node_map_it_end = node_to_node_map.end();
-      for ( ; node_map_it != node_map_it_end; ++node_map_it)
-        {
-          node_map_it->second += node_delta;
-        }
-      std::map<dof_id_type, std::vector<dof_id_type>>::iterator elem_map_it     = node_to_elems_map.begin();
-      std::map<dof_id_type, std::vector<dof_id_type>>::iterator elem_map_it_end = node_to_elems_map.end();
-      for ( ; elem_map_it != elem_map_it_end; ++elem_map_it)
-        {
-          std::size_t n_elems = elem_map_it->second.size();
-          for (std::size_t i=0; i<n_elems; i++)
-            {
-              (elem_map_it->second)[i] += elem_delta;
-            }
-        }
+      for (auto & pr : node_to_node_map)
+        pr.second += node_delta;
+
+      for (auto & pr : node_to_elems_map)
+        for (auto & entry : pr.second)
+          entry += elem_delta;
 
       // Copy mesh data. If we skip the call to find_neighbors(), the lists
       // of neighbors will be copied verbatim from the other mesh
       this->copy_nodes_and_elements(*other_mesh, skip_find_neighbors);
 
       // Decrement node IDs of mesh to return to original state
-      node_it  = other_mesh->nodes_begin();
-      node_end = other_mesh->nodes_end();
-      for (; node_it != node_end; ++node_it)
+      for (auto & nd : other_mesh->node_ptr_range())
         {
-          Node * nd = *node_it;
           dof_id_type new_id = nd->id() - node_delta;
           nd->set_id(new_id);
         }
@@ -1236,12 +1187,8 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
       // Container to catch boundary IDs passed back from BoundaryInfo.
       std::vector<boundary_id_type> bc_ids;
 
-      elem_it  = other_mesh->elements_begin();
-      elem_end = other_mesh->elements_end();
-      for (; elem_it != elem_end; ++elem_it)
+      for (auto & other_elem : other_mesh->element_ptr_range())
         {
-          Elem * other_elem = *elem_it;
-
           // Decrement elem IDs of other_mesh to return it to original state
           dof_id_type new_id = other_elem->id() - elem_delta;
           other_elem->set_id(new_id);
@@ -1314,21 +1261,19 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
   // Container to catch boundary IDs passed back from BoundaryInfo.
   std::vector<boundary_id_type> bc_ids;
 
-  std::map<dof_id_type, std::vector<dof_id_type>>::iterator elem_map_it     = node_to_elems_map.begin();
-  std::map<dof_id_type, std::vector<dof_id_type>>::iterator elem_map_it_end = node_to_elems_map.end();
   {
     LOG_SCOPE("stitch_meshes node updates", "ReplicatedMesh");
 
-    for ( ; elem_map_it != elem_map_it_end; ++elem_map_it)
+    for (const auto & pr : node_to_elems_map)
       {
-        dof_id_type target_node_id = elem_map_it->first;
+        dof_id_type target_node_id = pr.first;
         dof_id_type other_node_id = node_to_node_map[target_node_id];
         Node & target_node = this->node_ref(target_node_id);
 
-        std::size_t n_elems = elem_map_it->second.size();
+        std::size_t n_elems = pr.second.size();
         for (std::size_t i=0; i<n_elems; i++)
           {
-            dof_id_type elem_id = elem_map_it->second[i];
+            dof_id_type elem_id = pr.second[i];
             Elem * el = this->elem_ptr(elem_id);
 
             // find the local node index that we want to update
@@ -1343,18 +1288,16 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
       }
   }
 
-  std::map<dof_id_type, dof_id_type>::iterator node_map_it     = node_to_node_map.begin();
-  std::map<dof_id_type, dof_id_type>::iterator node_map_it_end = node_to_node_map.end();
   {
     LOG_SCOPE("stitch_meshes node deletion", "ReplicatedMesh");
-    for ( ; node_map_it != node_map_it_end; ++node_map_it)
+    for (const auto & pr : node_to_node_map)
       {
         // In the case that this==other_mesh, the two nodes might be the same (e.g. if
         // we're stitching a "sliver"), hence we need to skip node deletion in that case.
-        if ((this == other_mesh) && (node_map_it->second == node_map_it->first))
+        if ((this == other_mesh) && (pr.second == pr.first))
           continue;
 
-        dof_id_type this_node_id = node_map_it->second;
+        dof_id_type this_node_id = pr.second;
         this->delete_node( this->node_ptr(this_node_id) );
       }
   }
@@ -1370,15 +1313,13 @@ void ReplicatedMesh::stitching_helper (ReplicatedMesh * other_mesh,
     {
       LOG_SCOPE("stitch_meshes neighbor fixes", "ReplicatedMesh");
 
-      elem_map_it     = node_to_elems_map.begin();
-      elem_map_it_end = node_to_elems_map.end();
       std::set<dof_id_type> fixed_elems;
-      for ( ; elem_map_it != elem_map_it_end; ++elem_map_it)
+      for (const auto & pr : node_to_elems_map)
         {
-          std::size_t n_elems = elem_map_it->second.size();
+          std::size_t n_elems = pr.second.size();
           for (std::size_t i=0; i<n_elems; i++)
             {
-              dof_id_type elem_id = elem_map_it->second[i];
+              dof_id_type elem_id = pr.second[i];
               if (fixed_elems.find(elem_id) == fixed_elems.end())
                 {
                   Elem * el = this->elem_ptr(elem_id);
diff --git a/src/mesh/tecplot_io.C b/src/mesh/tecplot_io.C
index 16899f51fa4..7beab5e38b3 100644
--- a/src/mesh/tecplot_io.C
+++ b/src/mesh/tecplot_io.C
@@ -499,13 +499,12 @@ void TecplotIO::write_binary (const std::string & fname,
 
   // A zone for each subdomain
   bool firstzone=true;
-  for (std::set<subdomain_id_type>::const_iterator sbd_it=_subdomain_ids.begin();
-       sbd_it!=_subdomain_ids.end(); ++sbd_it)
+  for (const auto & sbd_id : _subdomain_ids)
     {
       // Copy the connectivity for this subdomain
       {
-        MeshBase::const_element_iterator       it  = the_mesh.active_subdomain_elements_begin (*sbd_it);
-        const MeshBase::const_element_iterator end = the_mesh.active_subdomain_elements_end   (*sbd_it);
+        MeshBase::const_element_iterator       it  = the_mesh.active_subdomain_elements_begin (sbd_id);
+        const MeshBase::const_element_iterator end = the_mesh.active_subdomain_elements_end   (sbd_id);
 
         unsigned int n_subcells_in_subdomain=0;
 
@@ -517,7 +516,7 @@ void TecplotIO::write_binary (const std::string & fname,
 
         unsigned int te = 0;
 
-        for (it  = the_mesh.active_subdomain_elements_begin (*sbd_it);
+        for (it  = the_mesh.active_subdomain_elements_begin (sbd_id);
              it != end; ++it)
           {
             std::vector<dof_id_type> conn;
@@ -543,7 +542,7 @@ void TecplotIO::write_binary (const std::string & fname,
           i_cell_max  = 0,
           j_cell_max  = 0,
           k_cell_max  = 0,
-          strand_id   = std::max(*sbd_it,static_cast<subdomain_id_type>(1)) + this->strand_offset(),
+          strand_id   = std::max(sbd_id, static_cast<subdomain_id_type>(1)) + this->strand_offset(),
           parent_zone = 0,
           is_block    = 1,
           num_face_connect   = 0,
@@ -559,7 +558,7 @@ void TecplotIO::write_binary (const std::string & fname,
         // zones will share from this one.
 
         // get the subdomain name from libMesh, if there is one.
-        std::string subdomain_name = the_mesh.subdomain_name(*sbd_it);
+        std::string subdomain_name = the_mesh.subdomain_name(sbd_id);
         std::ostringstream zone_name;
         zone_name << this->zone_title();
 
@@ -575,7 +574,7 @@ void TecplotIO::write_binary (const std::string & fname,
         else if (_subdomain_ids.size() > 1)
           {
             zone_name << "_";
-            zone_name << *sbd_it;
+            zone_name << sbd_id;
           }
 
         ierr = TECZNE112 (const_cast<char *>(zone_name.str().c_str()),
diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C
index c70f055fedf..b51a2a567c8 100644
--- a/src/mesh/unstructured_mesh.C
+++ b/src/mesh/unstructured_mesh.C
@@ -242,18 +242,12 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
 
   LOG_SCOPE("find_neighbors()", "Mesh");
 
-  const element_iterator el_end = this->elements_end();
-
   //TODO:[BSK] This should be removed later?!
   if (reset_current_list)
-    for (element_iterator el = this->elements_begin(); el != el_end; ++el)
-      {
-        Elem * e = *el;
-        for (auto s : e->side_index_range())
-          if (e->neighbor_ptr(s) != remote_elem ||
-              reset_remote_elements)
-            e->set_neighbor(s, libmesh_nullptr);
-      }
+    for (const auto & e : this->element_ptr_range())
+      for (auto s : e->side_index_range())
+        if (e->neighbor_ptr(s) != remote_elem || reset_remote_elements)
+          e->set_neighbor(s, libmesh_nullptr);
 
   // Find neighboring elements by first finding elements
   // with identical side keys and then check to see if they
@@ -271,10 +265,8 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
 
 
 
-    for (element_iterator el = this->elements_begin(); el != el_end; ++el)
+    for (const auto & element : this->element_ptr_range())
       {
-        Elem * element = *el;
-
         for (auto ms : element->side_index_range())
           {
           next_side:
diff --git a/src/mesh/xdr_io.C b/src/mesh/xdr_io.C
index 91f7fd7ed4a..856bed4750d 100644
--- a/src/mesh/xdr_io.C
+++ b/src/mesh/xdr_io.C
@@ -330,16 +330,13 @@ void XdrIO::write_serialized_subdomain_names(Xdr & io) const
       // return writable references in mesh_base, it's possible for the user to leave some entity names
       // blank.  We can't write those to the XDA file.
       new_header_id_type n_subdomain_names = 0;
-      std::map<subdomain_id_type, std::string>::const_iterator it_end = subdomain_map.end();
-      for (std::map<subdomain_id_type, std::string>::const_iterator it = subdomain_map.begin(); it != it_end; ++it)
-        {
-          if (!it->second.empty())
-            {
-              n_subdomain_names++;
-              subdomain_ids.push_back(it->first);
-              subdomain_names.push_back(it->second);
-            }
-        }
+      for (const auto & pr : subdomain_map)
+        if (!pr.second.empty())
+          {
+            n_subdomain_names++;
+            subdomain_ids.push_back(pr.first);
+            subdomain_names.push_back(pr.second);
+          }
 
       io.data(n_subdomain_names, "# subdomain id to name map");
       // Write out the ids and names in two vectors
@@ -1082,16 +1079,13 @@ void XdrIO::write_serialized_bcs_helper (Xdr & io, const new_header_id_type n_bc
           for (auto s : elem->side_index_range())
             {
               boundary_info.boundary_ids (elem, s, bc_ids);
-              for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
-                {
-                  const boundary_id_type bc_id = *id_it;
-                  if (bc_id != BoundaryInfo::invalid_id)
-                    {
-                      xfer_bcs.push_back (n_local_level_0_elem);
-                      xfer_bcs.push_back (s) ;
-                      xfer_bcs.push_back (bc_id);
-                    }
-                }
+              for (const auto & bc_id : bc_ids)
+                if (bc_id != BoundaryInfo::invalid_id)
+                  {
+                    xfer_bcs.push_back (n_local_level_0_elem);
+                    xfer_bcs.push_back (s) ;
+                    xfer_bcs.push_back (bc_id);
+                  }
             }
         }
       else if (bc_type == "edge")
@@ -1099,16 +1093,13 @@ void XdrIO::write_serialized_bcs_helper (Xdr & io, const new_header_id_type n_bc
           for (auto e : elem->edge_index_range())
             {
               boundary_info.edge_boundary_ids (elem, e, bc_ids);
-              for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
-                {
-                  const boundary_id_type bc_id = *id_it;
-                  if (bc_id != BoundaryInfo::invalid_id)
-                    {
-                      xfer_bcs.push_back (n_local_level_0_elem);
-                      xfer_bcs.push_back (e) ;
-                      xfer_bcs.push_back (bc_id);
-                    }
-                }
+              for (const auto & bc_id : bc_ids)
+                if (bc_id != BoundaryInfo::invalid_id)
+                  {
+                    xfer_bcs.push_back (n_local_level_0_elem);
+                    xfer_bcs.push_back (e) ;
+                    xfer_bcs.push_back (bc_id);
+                  }
             }
         }
       else if (bc_type == "shellface")
@@ -1116,16 +1107,13 @@ void XdrIO::write_serialized_bcs_helper (Xdr & io, const new_header_id_type n_bc
           for (unsigned short sf=0; sf<2; sf++)
             {
               boundary_info.shellface_boundary_ids (elem, sf, bc_ids);
-              for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
-                {
-                  const boundary_id_type bc_id = *id_it;
-                  if (bc_id != BoundaryInfo::invalid_id)
-                    {
-                      xfer_bcs.push_back (n_local_level_0_elem);
-                      xfer_bcs.push_back (sf) ;
-                      xfer_bcs.push_back (bc_id);
-                    }
-                }
+              for (const auto & bc_id : bc_ids)
+                if (bc_id != BoundaryInfo::invalid_id)
+                  {
+                    xfer_bcs.push_back (n_local_level_0_elem);
+                    xfer_bcs.push_back (sf) ;
+                    xfer_bcs.push_back (bc_id);
+                  }
             }
         }
       else
@@ -1221,15 +1209,12 @@ void XdrIO::write_serialized_nodesets (Xdr & io, const new_header_id_type n_node
   for (const auto & node : mesh.local_node_ptr_range())
     {
       boundary_info.boundary_ids (node, nodeset_ids);
-      for (std::vector<boundary_id_type>::const_iterator id_it=nodeset_ids.begin(); id_it!=nodeset_ids.end(); ++id_it)
-        {
-          const boundary_id_type bc_id = *id_it;
-          if (bc_id != BoundaryInfo::invalid_id)
-            {
-              xfer_bcs.push_back (node->id());
-              xfer_bcs.push_back (bc_id);
-            }
-        }
+      for (const auto & bc_id : nodeset_ids)
+        if (bc_id != BoundaryInfo::invalid_id)
+          {
+            xfer_bcs.push_back (node->id());
+            xfer_bcs.push_back (bc_id);
+          }
     }
 
   xfer_bcs.push_back(n_node);
@@ -1285,16 +1270,13 @@ void XdrIO::write_serialized_bc_names (Xdr & io, const BoundaryInfo & info, bool
       // return writable references in boundary_info, it's possible for the user to leave some entity names
       // blank.  We can't write those to the XDA file.
       new_header_id_type n_boundary_names = 0;
-      std::map<boundary_id_type, std::string>::const_iterator it_end = boundary_map.end();
-      for (std::map<boundary_id_type, std::string>::const_iterator it = boundary_map.begin(); it != it_end; ++it)
-        {
-          if (!it->second.empty())
-            {
-              n_boundary_names++;
-              boundary_ids.push_back(it->first);
-              boundary_names.push_back(it->second);
-            }
-        }
+      for (const auto & pr : boundary_map)
+        if (!pr.second.empty())
+          {
+            n_boundary_names++;
+            boundary_ids.push_back(pr.first);
+            boundary_names.push_back(pr.second);
+          }
 
       if (is_sideset)
         io.data(n_boundary_names, "# sideset id to name map");
diff --git a/src/parallel/parallel_elem.C b/src/parallel/parallel_elem.C
index 8a5216d898f..8b3cc77996e 100644
--- a/src/parallel/parallel_elem.C
+++ b/src/parallel/parallel_elem.C
@@ -295,9 +295,8 @@ Packing<const Elem *>::pack (const Elem * const & elem,
 
           *data_out++ =(bcs.size());
 
-          for (std::vector<boundary_id_type>::iterator bc_it=bcs.begin();
-               bc_it != bcs.end(); ++bc_it)
-            *data_out++ =(*bc_it);
+          for (const auto & bid : bcs)
+            *data_out++ = bid;
         }
 
       for (auto e : elem->edge_index_range())
@@ -306,9 +305,8 @@ Packing<const Elem *>::pack (const Elem * const & elem,
 
           *data_out++ =(bcs.size());
 
-          for (std::vector<boundary_id_type>::iterator bc_it=bcs.begin();
-               bc_it != bcs.end(); ++bc_it)
-            *data_out++ =(*bc_it);
+          for (const auto & bid : bcs)
+            *data_out++ = bid;
         }
 
       for (unsigned short sf=0; sf != 2; ++sf)
@@ -317,9 +315,8 @@ Packing<const Elem *>::pack (const Elem * const & elem,
 
           *data_out++ =(bcs.size());
 
-          for (std::vector<boundary_id_type>::iterator bc_it=bcs.begin();
-               bc_it != bcs.end(); ++bc_it)
-            *data_out++ =(*bc_it);
+          for (const auto & bid : bcs)
+            *data_out++ = bid;
         }
     }
 }
diff --git a/src/parallel/parallel_node.C b/src/parallel/parallel_node.C
index 1839189e851..d51cd7d7790 100644
--- a/src/parallel/parallel_node.C
+++ b/src/parallel/parallel_node.C
@@ -177,9 +177,8 @@ Packing<const Node *>::pack (const Node * const & node,
 
   *data_out++ =(bcs.size());
 
-  for (std::vector<boundary_id_type>::iterator bc_it=bcs.begin();
-       bc_it != bcs.end(); ++bc_it)
-    *data_out++ =(*bc_it);
+  for (const auto & bid : bcs)
+    *data_out++ = bid;
 }
 
 
diff --git a/src/quadrature/quadrature_composite.C b/src/quadrature/quadrature_composite.C
index 1a4739a1227..d5ccdbdc412 100644
--- a/src/quadrature/quadrature_composite.C
+++ b/src/quadrature/quadrature_composite.C
@@ -123,15 +123,14 @@ void QComposite<QSubCell>::add_subelem_values (const std::vector<Elem const *> &
   const std::vector<Real>  & subelem_weights = _lagrange_fe->get_JxW();
   const std::vector<Point> & subelem_points  = _lagrange_fe->get_xyz();
 
-  for (std::vector<Elem const *>::const_iterator it = subelem.begin();
-       it!=subelem.end(); ++it)
+  for (const auto & elem : subelem)
     {
       // tetgen seems to create 0-volume cells on occasion, but we *should*
       // be catching that appropriately now inside the ElemCutter class.
       // Just in case trap here, describe the error, and abort.
       libmesh_try
         {
-          _lagrange_fe->reinit(*it);
+          _lagrange_fe->reinit(elem);
           _weights.insert(_weights.end(),
                           subelem_weights.begin(), subelem_weights.end());
 
@@ -142,8 +141,8 @@ void QComposite<QSubCell>::add_subelem_values (const std::vector<Elem const *> &
         {
           libMesh::err << "ERROR: found a bad cut cell!\n";
 
-          for (unsigned int n=0; n<(*it)->n_nodes(); n++)
-            libMesh::err << (*it)->point(n) << std::endl;
+          for (unsigned int n=0; n<elem->n_nodes(); n++)
+            libMesh::err << elem->point(n) << std::endl;
 
           libmesh_error_msg("Tetgen may have created a 0-volume cell during Cutcell integration.");
         }
diff --git a/src/reduced_basis/rb_construction.C b/src/reduced_basis/rb_construction.C
index 93f9b272c27..3d481da8ec0 100644
--- a/src/reduced_basis/rb_construction.C
+++ b/src/reduced_basis/rb_construction.C
@@ -245,17 +245,11 @@ void RBConstruction::process_parameters_file (const std::string & parameters_fil
     }
 
   std::map<std::string,bool> log_scaling_in;
-  RBParameters::const_iterator it     = mu_min_in.begin();
-  RBParameters::const_iterator it_end = mu_min_in.end();
-  for ( ; it != it_end; ++it)
-    {
-      std::string param_name = it->first;
-
-      // For now, just set all entries to false.
-      // TODO: Implement a decent way to specify log-scaling true/false
-      // in the input text file
-      log_scaling_in[param_name] = false;
-    }
+  // For now, just set all entries to false.
+  // TODO: Implement a decent way to specify log-scaling true/false
+  // in the input text file
+  for (const auto & pr : mu_min_in)
+    log_scaling_in[pr.first] = false;
 
   // Set the parameters that have been read in
   set_rb_construction_parameters(n_training_samples,
@@ -334,18 +328,14 @@ void RBConstruction::print_info()
       libMesh::out << "RBThetaExpansion member is not set yet" << std::endl;
     }
   libMesh::out << "Number of parameters: " << get_n_params() << std::endl;
-  RBParameters::const_iterator it     = get_parameters().begin();
-  RBParameters::const_iterator it_end = get_parameters().end();
-  for ( ; it != it_end; ++it)
-    {
-      std::string param_name = it->first;
-      if (!is_discrete_parameter(param_name))
-        {
-          libMesh::out <<   "Parameter " << param_name
-                       << ": Min = " << get_parameter_min(param_name)
-                       << ", Max = " << get_parameter_max(param_name) << std::endl;
-        }
-    }
+  for (const auto & pr : get_parameters())
+    if (!is_discrete_parameter(pr.first))
+      {
+        libMesh::out <<   "Parameter " << pr.first
+                     << ": Min = " << get_parameter_min(pr.first)
+                     << ", Max = " << get_parameter_max(pr.first) << std::endl;
+      }
+
   print_discrete_parameter_values();
   libMesh::out << "n_training_samples: " << get_n_training_samples() << std::endl;
   libMesh::out << "quiet mode? " << is_quiet() << std::endl;
@@ -645,16 +635,8 @@ void RBConstruction::add_scaled_matrix_and_vector(Number scalar,
                   std::map<numeric_index_type, Number> rhs_values;
                   elem_assembly->get_nodal_rhs_values(rhs_values, *this, node);
 
-                  std::map<numeric_index_type, Number>::const_iterator it =
-                    rhs_values.begin();
-                  const std::map<numeric_index_type, Number>::const_iterator it_end =
-                    rhs_values.end();
-                  for ( ; it != it_end; ++it)
-                    {
-                      numeric_index_type dof_index = it->first;
-                      Number value = it->second;
-                      input_vector->add( dof_index, value);
-                    }
+                  for (const auto & pr : rhs_values)
+                    input_vector->add(pr.first, pr.second);
                 }
 #endif
             }
@@ -788,12 +770,8 @@ void RBConstruction::add_scaled_matrix_and_vector(Number scalar,
               for (unsigned int var1=0; var1<n_vars(); var1++)
                 {
                   ConstCouplingRow ccr(var1, *coupling_matrix);
-                  ConstCouplingRow::const_iterator end = ccr.end();
-                  for (ConstCouplingRow::const_iterator it =
-                         ccr.begin(); it != end; ++it)
+                  for (const auto & var2 : ccr)
                     {
-                      unsigned int var2 = *it;
-
                       unsigned int sub_m = context.get_elem_jacobian( var1, var2 ).m();
                       unsigned int sub_n = context.get_elem_jacobian( var1, var2 ).n();
                       DenseMatrix<Number> sub_jac(sub_m, sub_n);
diff --git a/src/solution_transfer/dtk_adapter.C b/src/solution_transfer/dtk_adapter.C
index 719950437c0..891f870550f 100644
--- a/src/solution_transfer/dtk_adapter.C
+++ b/src/solution_transfer/dtk_adapter.C
@@ -57,11 +57,9 @@ DTKAdapter::DTKAdapter(Teuchos::RCP<const Teuchos::Comm<int>> in_comm, EquationS
   {
     unsigned int i = 0;
 
-    for (std::set<unsigned int>::iterator it = semi_local_nodes.begin();
-         it != semi_local_nodes.end();
-         ++it)
+    for (const auto & id : semi_local_nodes)
       {
-        const Node & node = mesh.node_ref(*it);
+        const Node & node = mesh.node_ref(id);
 
         vertices[i] = node.id();
 
@@ -211,7 +209,7 @@ DTKAdapter::update_variable_values(std::string var_name)
 
   unsigned int i=0;
   // Loop over the values (one for each node) and assign the value of this variable at each node
-  for (FieldContainerType::iterator it=values->begin(); it != values->end(); ++it)
+  for (const auto & value : *values)
     {
       unsigned int node_num = vertices[i];
       const Node & node = mesh.node_ref(node_num);
@@ -220,7 +218,7 @@ DTKAdapter::update_variable_values(std::string var_name)
         {
           // The 0 is for the component... this only works for LAGRANGE!
           dof_id_type dof = node.dof_number(sys->number(), var_num, 0);
-          sys->solution->set(dof, *it);
+          sys->solution->set(dof, value);
         }
 
       i++;
diff --git a/src/solution_transfer/dtk_solution_transfer.C b/src/solution_transfer/dtk_solution_transfer.C
index bac524fea1b..45795ed5cf5 100644
--- a/src/solution_transfer/dtk_solution_transfer.C
+++ b/src/solution_transfer/dtk_solution_transfer.C
@@ -56,15 +56,11 @@ DTKSolutionTransfer::DTKSolutionTransfer(const libMesh::Parallel::Communicator &
 
 DTKSolutionTransfer::~DTKSolutionTransfer()
 {
-  for (std::map<EquationSystems *, DTKAdapter *>::iterator it = adapters.begin();
-       it != adapters.end();
-       ++it)
-    delete it->second;
-
-  for (std::map<std::pair<EquationSystems *, EquationSystems *>, shared_domain_map_type *>::iterator it = dtk_maps.begin();
-       it != dtk_maps.end();
-       ++it)
-    delete it->second;
+  for (auto & pr : adapters)
+    delete pr.second;
+
+  for (auto & pr : dtk_maps)
+    delete pr.second;
 }
 
 void
diff --git a/src/solution_transfer/meshfree_interpolation.C b/src/solution_transfer/meshfree_interpolation.C
index 33345f17621..0510f2946ac 100644
--- a/src/solution_transfer/meshfree_interpolation.C
+++ b/src/solution_transfer/meshfree_interpolation.C
@@ -223,10 +223,8 @@ void InverseDistanceInterpolation<KDDim>::interpolate_field_data (const std::vec
     std::vector<size_t> ret_index(num_results);
     std::vector<Real>   ret_dist_sqr(num_results);
 
-    for (std::vector<Point>::const_iterator tgt_it=tgt_pts.begin();
-         tgt_it != tgt_pts.end(); ++tgt_it)
+    for (const auto & tgt : tgt_pts)
       {
-        const Point & tgt(*tgt_it);
         const Real query_pt[] = { tgt(0), tgt(1), tgt(2) };
 
         _kd_tree->knnSearch(&query_pt[0], num_results, &ret_index[0], &ret_dist_sqr[0]);
diff --git a/src/solvers/petsc_auto_fieldsplit.C b/src/solvers/petsc_auto_fieldsplit.C
index 82af5470519..41974dc56bc 100644
--- a/src/solvers/petsc_auto_fieldsplit.C
+++ b/src/solvers/petsc_auto_fieldsplit.C
@@ -101,11 +101,8 @@ void petsc_auto_fieldsplit (PC my_pc,
         }
     }
 
-  for (std::map<std::string, std::vector<dof_id_type>>::const_iterator
-         i = group_indices.begin(); i != group_indices.end(); ++i)
-    {
-      indices_to_fieldsplit(sys.comm(), i->second, my_pc, i->first);
-    }
+  for (const auto & pr : group_indices)
+    indices_to_fieldsplit(sys.comm(), pr.second, my_pc, pr.first);
 }
 
 } // namespace libMesh
diff --git a/src/solvers/petscdmlibmeshimpl.C b/src/solvers/petscdmlibmeshimpl.C
index e002e1d6b57..d7f89cd7670 100644
--- a/src/solvers/petscdmlibmeshimpl.C
+++ b/src/solvers/petscdmlibmeshimpl.C
@@ -179,10 +179,11 @@ PetscErrorCode DMlibMeshGetBlocks(DM dm, PetscInt * n, char *** blocknames)
   if (!blocknames) PetscFunctionReturn(0);
   ierr = PetscMalloc(*n*sizeof(char *), blocknames); CHKERRQ(ierr);
   i = 0;
-  for (std::map<std::string, unsigned int>::const_iterator it = dlm->blockids->begin(); it != dlm->blockids->end(); ++it){
-    ierr = PetscStrallocpy(it->first.c_str(), *blocknames+i); CHKERRQ(ierr);
-    ++i;
-  }
+  for (const auto & pr : *(dlm->blockids))
+    {
+      ierr = PetscStrallocpy(pr.first.c_str(), *blocknames+i); CHKERRQ(ierr);
+      ++i;
+    }
   PetscFunctionReturn(0);
 }
 
@@ -203,10 +204,11 @@ PetscErrorCode DMlibMeshGetVariables(DM dm, PetscInt * n, char *** varnames)
   if (!varnames) PetscFunctionReturn(0);
   ierr = PetscMalloc(*n*sizeof(char *), varnames); CHKERRQ(ierr);
   i = 0;
-  for (std::map<std::string, unsigned int>::const_iterator it = dlm->varids->begin(); it != dlm->varids->end(); ++it){
-    ierr = PetscStrallocpy(it->first.c_str(), *varnames+i); CHKERRQ(ierr);
-    ++i;
-  }
+  for (const auto & pr : *(dlm->varids))
+    {
+      ierr = PetscStrallocpy(pr.first.c_str(), *varnames+i); CHKERRQ(ierr);
+      ++i;
+    }
   PetscFunctionReturn(0);
 }
 
@@ -282,8 +284,7 @@ static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len,
     std::map<std::string, unsigned int>  dvarids;
     std::map<unsigned int, std::string>  dvarnames;
     unsigned int                         dvcount = 0;
-    for (std::set<unsigned int>::const_iterator dvit = (*dlm->decomposition)[d].begin(); dvit != (*dlm->decomposition)[d].end(); ++dvit){
-      unsigned int v = *dvit;
+    for (const auto & v : (*dlm->decomposition)[d]) {
       std::string vname = (*dlm->varnames)[v];
       dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
       dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
@@ -291,11 +292,11 @@ static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len,
       else   dname += "_" + vname;
       ++dvcount;
       if (!islist) continue;
-      /* Iterate only over this DM's blocks. */
-      for (std::map<std::string, unsigned int>::const_iterator bit = dlm->blockids->begin(); bit != dlm->blockids->end(); ++bit) {
-        unsigned int b = bit->second;
-        MeshBase::const_element_iterator el     = sys->get_mesh().active_local_subdomain_elements_begin(b);
-        MeshBase::const_element_iterator end_el = sys->get_mesh().active_local_subdomain_elements_end(b);
+      // Iterate only over this DM's blocks.
+      for (const auto & pr : *(dlm->blockids)) {
+        MeshBase::const_element_iterator
+          el     = sys->get_mesh().active_local_subdomain_elements_begin(pr.second),
+          end_el = sys->get_mesh().active_local_subdomain_elements_end(pr.second);
         for ( ; el != end_el; ++el) {
           const Elem * elem = *el;
           //unsigned int e_subdomain = elem->subdomain_id();
@@ -304,7 +305,7 @@ static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len,
           dofmap.dof_indices(elem, evindices, v);
           for (unsigned int i = 0; i < evindices.size(); ++i) {
             numeric_index_type dof = evindices[i];
-            if (dof >= dofmap.first_dof() && dof < dofmap.end_dof()) /* might want to use variable_first/last_local_dof instead */
+            if (dof >= dofmap.first_dof() && dof < dofmap.end_dof()) // might want to use variable_first/last_local_dof instead
               dindices.insert(dof);
           }
         }
@@ -318,8 +319,8 @@ static PetscErrorCode  DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len,
       PetscInt * darray;
       ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
       numeric_index_type i = 0;
-      for (std::set<numeric_index_type>::const_iterator it = dindices.begin(); it != dindices.end(); ++it) {
-        darray[i] = *it;
+      for (const auto & id : dindices) {
+        darray[i] = id;
         ++i;
       }
       ierr = ISCreateGeneral(((PetscObject)dm)->comm, dindices.size(),darray, PETSC_OWN_POINTER, &dis); CHKERRQ(ierr);
@@ -387,8 +388,7 @@ static PetscErrorCode  DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len
     std::map<std::string, unsigned int>  dblockids;
     std::map<unsigned int,std::string>   dblocknames;
     unsigned int                         dbcount = 0;
-    for (std::set<unsigned int>::const_iterator bit = (*dlm->decomposition)[d].begin(); bit != (*dlm->decomposition)[d].end(); ++bit){
-      unsigned int b = *bit;
+    for (const auto & b : (*dlm->decomposition)[d]) {
       std::string bname = (*dlm->blocknames)[b];
       dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
       dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
@@ -401,14 +401,12 @@ static PetscErrorCode  DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len
       for ( ; el != end_el; ++el) {
         const Elem * elem = *el;
         std::vector<numeric_index_type> evindices;
-        /* Iterate only over this DM's variables. */
-        for (std::map<std::string, unsigned int>::const_iterator vit = dlm->varids->begin(); vit != dlm->varids->end(); ++vit) {
-          unsigned int v = vit->second;
+        // Iterate only over this DM's variables.
+        for (const auto & pr : *(dlm->varids)) {
           // Get the degree of freedom indices for the given variable off the current element.
-          sys->get_dof_map().dof_indices(elem, evindices, v);
-          for (unsigned int i = 0; i < evindices.size(); ++i) {
-            numeric_index_type dof = evindices[i];
-            if (dof >= sys->get_dof_map().first_dof() && dof < sys->get_dof_map().end_dof()) /* might want to use variable_first/last_local_dof instead */
+          sys->get_dof_map().dof_indices(elem, evindices, pr.second);
+          for (const auto & dof : evindices) {
+            if (dof >= sys->get_dof_map().first_dof() && dof < sys->get_dof_map().end_dof()) // might want to use variable_first/last_local_dof instead
               dindices.insert(dof);
           }
         }
@@ -422,8 +420,8 @@ static PetscErrorCode  DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len
       IS dis;
       ierr = PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray); CHKERRQ(ierr);
       numeric_index_type i = 0;
-      for (std::set<numeric_index_type>::const_iterator it = dindices.begin(); it != dindices.end(); ++it) {
-        darray[i] = *it;
+      for (const auto & id : dindices) {
+        darray[i] = id;
         ++i;
       }
       ierr = ISCreateGeneral(((PetscObject)dm)->comm, dindices.size(),darray, PETSC_OWN_POINTER, &dis); CHKERRQ(ierr);
@@ -512,12 +510,12 @@ PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, Pets
       }
     }
   }
-  else { /* Empty splits indicate default: split all variables with one per split. */
+  else { // Empty splits indicate default: split all variables with one per split.
     PetscInt d = 0;
-    for (std::map<std::string, unsigned int>::const_iterator vit = ddlm->varids->begin(); vit != ddlm->varids->end(); ++vit) {
+    for (const auto & pr : (*ddlm->varids)) {
       ddlm->decomposition->push_back(std::set<unsigned int>());
-      unsigned int vid = vit->second;
-      std::string vname = vit->first;
+      unsigned int vid = pr.second;
+      std::string vname = pr.first;
       (*ddlm->decomposition)[d].insert(vid);
       ++d;
     }
@@ -565,12 +563,12 @@ PetscErrorCode DMlibMeshCreateDomainDecompositionDM(DM dm, PetscInt dnumber, Pet
       }
     }
   }
-  else { /* Empty splits indicate default: split all blocks with one per split. */
+  else { // Empty splits indicate default: split all blocks with one per split.
     PetscInt d = 0;
-    for (std::map<std::string, unsigned int>::const_iterator bit = ddlm->blockids->begin(); bit != ddlm->blockids->end(); ++bit) {
+    for (const auto & pr : (*ddlm->blockids)) {
       ddlm->decomposition->push_back(std::set<unsigned int>());
-      unsigned int bid = bit->second;
-      std::string bname = bit->first;
+      unsigned int bid = pr.second;
+      std::string bname = pr.first;
       (*ddlm->decomposition)[d].insert(bid);
       ++d;
     }
diff --git a/src/systems/dg_fem_context.C b/src/systems/dg_fem_context.C
index 887bc632b6c..37e0964c418 100644
--- a/src/systems/dg_fem_context.C
+++ b/src/systems/dg_fem_context.C
@@ -91,11 +91,9 @@ void DGFEMContext::neighbor_side_fe_reinit ()
   // the quadrature points on the current side
   std::vector<Point> qface_side_points;
   std::vector<Point> qface_neighbor_points;
-  std::map<FEType, std::unique_ptr<FEAbstract>>::iterator local_fe_end = _neighbor_side_fe.end();
-  for (std::map<FEType, std::unique_ptr<FEAbstract>>::iterator i = _neighbor_side_fe.begin();
-       i != local_fe_end; ++i)
+  for (auto & pr : _neighbor_side_fe)
     {
-      FEType neighbor_side_fe_type = i->first;
+      FEType neighbor_side_fe_type = pr.first;
       FEAbstract * side_fe = _side_fe[this->get_dim()][neighbor_side_fe_type].get();
       qface_side_points = side_fe->get_xyz();
 
@@ -105,7 +103,7 @@ void DGFEMContext::neighbor_side_fe_reinit ()
                                 qface_side_points,
                                 qface_neighbor_points);
 
-      i->second->reinit(&get_neighbor(), &qface_neighbor_points);
+      pr.second->reinit(&get_neighbor(), &qface_neighbor_points);
     }
 
   // Set boolean flag to indicate that the DG terms are active on this element
diff --git a/src/systems/diff_system.C b/src/systems/diff_system.C
index a4a94bc932e..7e9ea7bcf03 100644
--- a/src/systems/diff_system.C
+++ b/src/systems/diff_system.C
@@ -200,10 +200,9 @@ void DifferentiableSystem::add_second_order_dot_vars()
   const std::set<unsigned int> & second_order_vars = this->get_second_order_vars();
   if (!second_order_vars.empty())
     {
-      for (std::set<unsigned int>::const_iterator var_it = second_order_vars.begin();
-           var_it != second_order_vars.end(); ++var_it)
+      for (const auto & var_id : second_order_vars)
         {
-          const Variable & var = this->variable(*var_it);
+          const Variable & var = this->variable(var_id);
           std::string new_var_name = std::string("dot_")+var.name();
 
           unsigned int v_var_idx;
@@ -213,14 +212,14 @@ void DifferentiableSystem::add_second_order_dot_vars()
           else
             v_var_idx = this->add_variable( new_var_name, var.type(), &var.active_subdomains() );
 
-          _second_order_dot_vars.insert( std::pair<unsigned int,unsigned int>(*var_it,v_var_idx) );
+          _second_order_dot_vars.insert(std::pair<unsigned int, unsigned int>(var_id, v_var_idx));
 
           // The new velocities are time evolving variables of first order
           this->time_evolving( v_var_idx, 1 );
 
           // And if there are any boundary conditions set on the second order
           // variable, we also need to set it on its velocity variable.
-          this->add_dot_var_dirichlet_bcs( *var_it, v_var_idx );
+          this->add_dot_var_dirichlet_bcs(var_id, v_var_idx);
         }
     }
 }
@@ -332,15 +331,9 @@ bool DifferentiableSystem::have_first_order_scalar_vars() const
   bool have_first_order_scalar_vars = false;
 
   if (this->have_first_order_vars())
-    {
-      for (std::set<unsigned int>::const_iterator var_it = this->get_first_order_vars().begin();
-           var_it != this->get_first_order_vars().end();
-           ++var_it)
-        {
-          if (this->variable(*var_it).type().family == SCALAR)
-            have_first_order_scalar_vars = true;
-        }
-    }
+    for (const auto & var : this->get_first_order_vars())
+      if (this->variable(var).type().family == SCALAR)
+        have_first_order_scalar_vars = true;
 
   return have_first_order_scalar_vars;
 }
@@ -350,15 +343,9 @@ bool DifferentiableSystem::have_second_order_scalar_vars() const
   bool have_second_order_scalar_vars = false;
 
   if (this->have_second_order_vars())
-    {
-      for (std::set<unsigned int>::const_iterator var_it = this->get_second_order_vars().begin();
-           var_it != this->get_second_order_vars().end();
-           ++var_it)
-        {
-          if (this->variable(*var_it).type().family == SCALAR)
-            have_second_order_scalar_vars = true;
-        }
-    }
+    for (const auto & var : this->get_second_order_vars())
+      if (this->variable(var).type().family == SCALAR)
+        have_second_order_scalar_vars = true;
 
   return have_second_order_scalar_vars;
 }
diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C
index 64f0819006d..3534cacf5d0 100644
--- a/src/systems/equation_systems.C
+++ b/src/systems/equation_systems.C
@@ -1135,16 +1135,14 @@ std::string EquationSystems::get_info () const
   //       oss << "  n_parameters()=" << this->n_parameters() << '\n';
   //       oss << "   Parameters:\n";
 
-  //       for (std::map<std::string, Real>::const_iterator
-  //      param = _parameters.begin(); param != _parameters.end();
-  //    ++param)
-  // oss << "    "
-  //     << "\""
-  //     << param->first
-  //     << "\""
-  //     << "="
-  //     << param->second
-  //     << '\n';
+  //       for (const auto & pr : _parameters)
+  //         oss << "    "
+  //             << "\""
+  //             << pr.first
+  //             << "\""
+  //             << "="
+  //             << pr.second
+  //             << '\n';
   //     }
 
   return oss.str();
diff --git a/src/systems/equation_systems_io.C b/src/systems/equation_systems_io.C
index fc1a91b187c..10fc34ee767 100644
--- a/src/systems/equation_systems_io.C
+++ b/src/systems/equation_systems_io.C
@@ -335,22 +335,19 @@ void EquationSystems::_read_impl (const std::string & name,
 
       Xdr local_io (read_parallel_files ? local_file_name(this->processor_id(),name) : "", mode);
 
-      std::vector<std::pair<std::string, System *>>::iterator
-        pos = xda_systems.begin();
-
-      for (; pos != xda_systems.end(); ++pos)
+      for (auto & pr : xda_systems)
         if (read_legacy_format)
           {
             libmesh_deprecated();
 #ifdef LIBMESH_ENABLE_DEPRECATED
-            pos->second->read_legacy_data (io, read_additional_data);
+            pr.second->read_legacy_data (io, read_additional_data);
 #endif
           }
         else
           if (read_parallel_files)
-            pos->second->read_parallel_data<InValType>   (local_io, read_additional_data);
+            pr.second->read_parallel_data<InValType>   (local_io, read_additional_data);
           else
-            pos->second->read_serialized_data<InValType> (io, read_additional_data);
+            pr.second->read_serialized_data<InValType> (io, read_additional_data);
 
 
       // Undo the temporary numbering.
@@ -485,11 +482,9 @@ void EquationSystems::write(const std::string & name,
     const unsigned int proc_id = this->processor_id();
 
     unsigned int n_sys = 0;
-    for (std::map<std::string, System *>::const_iterator pos = _systems.begin();
-         pos != _systems.end(); ++pos)
-      {
-        if (! pos->second->hide_output()) n_sys++;
-      }
+    for (auto & pr : _systems)
+      if (!pr.second->hide_output())
+        n_sys++;
 
     // set the version number in the Xdr object
     io.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
@@ -517,17 +512,16 @@ void EquationSystems::write(const std::string & name,
         // Write the number of equation systems
         io.data (n_sys, "# No. of Equation Systems");
 
-        for (std::map<std::string, System *>::const_iterator pos = _systems.begin();
-             pos != _systems.end(); ++pos)
+        for (auto & pr : _systems)
           {
             // Ignore this system if it has been marked as hidden
-            if (pos->second->hide_output()) continue;
+            if (pr.second->hide_output()) continue;
 
             // 3.)
             // Write the name of the sys_num-th system
             {
-              const unsigned int sys_num = pos->second->number();
-              std::string sys_name       = pos->first;
+              const unsigned int sys_num = pr.second->number();
+              std::string sys_name       = pr.first;
 
               comment =  "# Name, System No. ";
               std::sprintf(buf, "%u", sys_num);
@@ -539,8 +533,8 @@ void EquationSystems::write(const std::string & name,
             // 4.)
             // Write the type of system handled
             {
-              const unsigned int sys_num = pos->second->number();
-              std::string sys_type       = pos->second->system_type();
+              const unsigned int sys_num = pr.second->number();
+              std::string sys_type       = pr.second->system_type();
 
               comment =  "# Type, System No. ";
               std::sprintf(buf, "%u", sys_num);
@@ -551,7 +545,7 @@ void EquationSystems::write(const std::string & name,
 
             // 5.) - 9.)
             // Let System::write_header() do the job
-            pos->second->write_header (io, version, write_additional_data);
+            pr.second->write_header (io, version, write_additional_data);
           }
       }
 
@@ -562,17 +556,16 @@ void EquationSystems::write(const std::string & name,
         // open a parallel buffer if warranted.
         Xdr local_io (write_parallel_files ? local_file_name(this->processor_id(),name) : "", mode);
 
-        for (std::map<std::string, System *>::const_iterator pos = _systems.begin();
-             pos != _systems.end(); ++pos)
+        for (auto & pr : _systems)
           {
             // Ignore this system if it has been marked as hidden
-            if (pos->second->hide_output()) continue;
+            if (pr.second->hide_output()) continue;
 
             // 10.) + 11.)
             if (write_parallel_files)
-              pos->second->write_parallel_data (local_io,write_additional_data);
+              pr.second->write_parallel_data (local_io,write_additional_data);
             else
-              pos->second->write_serialized_data (io,write_additional_data);
+              pr.second->write_serialized_data (io,write_additional_data);
           }
       }
   }
diff --git a/src/systems/fem_context.C b/src/systems/fem_context.C
index 229bb126cd2..31957cf912a 100644
--- a/src/systems/fem_context.C
+++ b/src/systems/fem_context.C
@@ -124,11 +124,8 @@ void FEMContext::init_internal_data(const System & sys)
     // SCALAR FEs have dimension 0 by assumption
     _elem_dims.insert(0);
 
-  for (std::set<unsigned char>::const_iterator dim_it = _elem_dims.begin();
-       dim_it != _elem_dims.end(); ++dim_it)
+  for (const auto & dim : _elem_dims)
     {
-      const unsigned char dim = *dim_it;
-
       // Create an adequate quadrature rule
       _element_qrule[dim] =
         hardest_fe_type.default_quadrature_rule(dim, _extra_quadrature_order);
@@ -1360,15 +1357,13 @@ void FEMContext::elem_fe_reinit(const std::vector<Point> * const pts)
 
   libmesh_assert( !_element_fe[dim].empty() );
 
-  std::map<FEType, std::unique_ptr<FEAbstract>>::iterator local_fe_end = _element_fe[dim].end();
-  for (std::map<FEType, std::unique_ptr<FEAbstract>>::iterator i = _element_fe[dim].begin();
-       i != local_fe_end; ++i)
+  for (const auto & pr : _element_fe[dim])
     {
       if (this->has_elem())
-        i->second->reinit(&(this->get_elem()), pts);
+        pr.second->reinit(&(this->get_elem()), pts);
       else
         // If !this->has_elem(), then we assume we are dealing with a SCALAR variable
-        i->second->reinit(libmesh_nullptr);
+        pr.second->reinit(libmesh_nullptr);
     }
 }
 
@@ -1383,12 +1378,8 @@ void FEMContext::side_fe_reinit ()
 
   libmesh_assert( !_side_fe[dim].empty() );
 
-  std::map<FEType, std::unique_ptr<FEAbstract>>::iterator local_fe_end = _side_fe[dim].end();
-  for (std::map<FEType, std::unique_ptr<FEAbstract>>::iterator i = _side_fe[dim].begin();
-       i != local_fe_end; ++i)
-    {
-      i->second->reinit(&(this->get_elem()), this->get_side());
-    }
+  for (auto & pr : _side_fe[dim])
+    pr.second->reinit(&(this->get_elem()), this->get_side());
 }
 
 
@@ -1399,12 +1390,8 @@ void FEMContext::edge_fe_reinit ()
 
   // Initialize all the interior FE objects on elem/edge.
   // Logging of FE::reinit is done in the FE functions
-  std::map<FEType, std::unique_ptr<FEAbstract>>::iterator local_fe_end = _edge_fe.end();
-  for (std::map<FEType, std::unique_ptr<FEAbstract>>::iterator i = _edge_fe.begin();
-       i != local_fe_end; ++i)
-    {
-      i->second->edge_reinit(&(this->get_elem()), this->get_edge());
-    }
+  for (auto & pr : _edge_fe)
+    pr.second->edge_reinit(&(this->get_elem()), this->get_edge());
 }
 
 
diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C
index cc12c58bee8..c7f5ecbb942 100644
--- a/src/systems/fem_system.C
+++ b/src/systems/fem_system.C
@@ -371,10 +371,9 @@ public:
     FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
     _sys.init_context(_femcontext);
 
-    for (ConstElemRange::const_iterator elem_it = range.begin();
-         elem_it != range.end(); ++elem_it)
+    for (const auto & elem : range)
       {
-        Elem * el = const_cast<Elem *>(*elem_it);
+        Elem * el = const_cast<Elem *>(elem);
 
         _femcontext.pre_fe_reinit(_sys, el);
         _femcontext.elem_fe_reinit();
@@ -413,10 +412,9 @@ public:
     FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
     _sys.init_context(_femcontext);
 
-    for (ConstElemRange::const_iterator elem_it = range.begin();
-         elem_it != range.end(); ++elem_it)
+    for (const auto & elem : range)
       {
-        Elem * el = const_cast<Elem *>(*elem_it);
+        Elem * el = const_cast<Elem *>(elem);
         _femcontext.pre_fe_reinit(_sys, el);
 
         // Optionally initialize all the interior FE objects on elem.
@@ -492,10 +490,9 @@ public:
     if (have_some_heterogenous_qoi_bc)
       _sys.init_context(_femcontext);
 
-    for (ConstElemRange::const_iterator elem_it = range.begin();
-         elem_it != range.end(); ++elem_it)
+    for (const auto & elem : range)
       {
-        Elem * el = const_cast<Elem *>(*elem_it);
+        Elem * el = const_cast<Elem *>(elem);
 
         _femcontext.pre_fe_reinit(_sys, el);
 
@@ -628,10 +625,9 @@ public:
     if (have_some_heterogenous_qoi_bc)
       _sys.init_context(_femcontext);
 
-    for (ConstElemRange::const_iterator elem_it = range.begin();
-         elem_it != range.end(); ++elem_it)
+    for (const auto & elem : range)
       {
-        Elem * el = const_cast<Elem *>(*elem_it);
+        Elem * el = const_cast<Elem *>(elem);
 
         _femcontext.pre_fe_reinit(_sys, el);
 
diff --git a/src/systems/implicit_system.C b/src/systems/implicit_system.C
index b53d2ab3a5d..5eae91d6f0c 100644
--- a/src/systems/implicit_system.C
+++ b/src/systems/implicit_system.C
@@ -70,12 +70,11 @@ void ImplicitSystem::clear ()
 
   // clear any user-added matrices
   {
-    for (matrices_iterator pos = _matrices.begin();
-         pos != _matrices.end(); ++pos)
+    for (auto & pr : _matrices)
       {
-        pos->second->clear ();
-        delete pos->second;
-        pos->second = libmesh_nullptr;
+        pr.second->clear ();
+        delete pr.second;
+        pr.second = libmesh_nullptr;
       }
 
     _matrices.clear();
@@ -94,9 +93,8 @@ void ImplicitSystem::init_data ()
   Parent::init_data();
 
   // Clear any existing matrices
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
-    pos->second->clear();
+  for (auto & pr : _matrices)
+    pr.second->clear();
 
   // Initialize the matrices for the system
   this->init_matrices ();
@@ -121,10 +119,9 @@ void ImplicitSystem::init_matrices ()
   _can_add_matrices = false;
 
   // Tell the matrices about the dof map, and vice versa
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
+  for (auto & pr : _matrices)
     {
-      SparseMatrix<Number> & m = *(pos->second);
+      SparseMatrix<Number> & m = *(pr.second);
       libmesh_assert (!m.initialized());
 
       // We want to allow repeated init() on systems, but we don't
@@ -139,14 +136,12 @@ void ImplicitSystem::init_matrices ()
   dof_map.compute_sparsity (this->get_mesh());
 
   // Initialize matrices
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
-    pos->second->init ();
+  for (auto & pr : _matrices)
+    pr.second->init ();
 
   // Set the additional matrices to 0.
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
-    pos->second->zero ();
+  for (auto & pr : _matrices)
+    pr.second->zero ();
 }
 
 
@@ -160,11 +155,10 @@ void ImplicitSystem::reinit ()
   DofMap & dof_map = this->get_dof_map();
 
   // Clear the matrices
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
+  for (auto & pr : _matrices)
     {
-      pos->second->clear();
-      pos->second->attach_dof_map (dof_map);
+      pr.second->clear();
+      pr.second->attach_dof_map (dof_map);
     }
 
   // Clear the sparsity pattern
@@ -176,14 +170,12 @@ void ImplicitSystem::reinit ()
   dof_map.compute_sparsity (this->get_mesh());
 
   // Initialize matrices
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
-    pos->second->init ();
+  for (auto & pr : _matrices)
+    pr.second->init ();
 
   // Set the additional matrices to 0.
-  for (matrices_iterator pos = _matrices.begin();
-       pos != _matrices.end(); ++pos)
-    pos->second->zero ();
+  for (auto & pr : _matrices)
+    pr.second->zero ();
 }
 
 
diff --git a/src/systems/qoi_set.C b/src/systems/qoi_set.C
index 05c092ddf0f..597d9468ad5 100644
--- a/src/systems/qoi_set.C
+++ b/src/systems/qoi_set.C
@@ -46,15 +46,13 @@ unsigned int QoISet::size (const System & sys) const
 void QoISet::add_indices(const std::vector<unsigned int> & indices)
 {
   unsigned int max_size = 0;
-  for (std::vector<unsigned int>::const_iterator i = indices.begin();
-       i != indices.end(); ++i)
-    max_size = std::max(max_size, *i + 1);
+  for (const auto & i : indices)
+    max_size = std::max(max_size, i + 1);
 
   _indices.resize(max_size);
 
-  for (std::vector<unsigned int>::const_iterator i = indices.begin();
-       i != indices.end(); ++i)
-    _indices[*i] = true;
+  for (const auto & i : indices)
+    _indices[i] = true;
 }
 
 
@@ -62,9 +60,8 @@ void QoISet::add_indices(const std::vector<unsigned int> & indices)
 inline
 void QoISet::remove_indices(const std::vector<unsigned int> & indices)
 {
-  for (std::vector<unsigned int>::const_iterator i = indices.begin();
-       i != indices.end(); ++i)
-    _indices[*i] = false;
+  for (const auto & i : indices)
+    _indices[i] = false;
 }
 
 } // namespace libMesh
diff --git a/src/systems/system.C b/src/systems/system.C
index 02203fcc2bf..5afed42fced 100644
--- a/src/systems/system.C
+++ b/src/systems/system.C
@@ -214,11 +214,11 @@ void System::clear ()
 
   // clear any user-added vectors
   {
-    for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
+    for (auto & pr : _vectors)
       {
-        pos->second->clear ();
-        delete pos->second;
-        pos->second = libmesh_nullptr;
+        pr.second->clear ();
+        delete pr.second;
+        pr.second = libmesh_nullptr;
       }
 
     _vectors.clear();
@@ -291,28 +291,28 @@ void System::init_data ()
   _is_initialized = true;
 
   // initialize & zero other vectors, if necessary
-  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
+  for (auto & pr : _vectors)
     {
-      ParallelType type = _vector_types[pos->first];
+      ParallelType type = _vector_types[pr.first];
 
       if (type == GHOSTED)
         {
 #ifdef LIBMESH_ENABLE_GHOSTED
-          pos->second->init (this->n_dofs(), this->n_local_dofs(),
-                             _dof_map->get_send_list(), false,
-                             GHOSTED);
+          pr.second->init (this->n_dofs(), this->n_local_dofs(),
+                           _dof_map->get_send_list(), false,
+                           GHOSTED);
 #else
           libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
 #endif
         }
       else if (type == SERIAL)
         {
-          pos->second->init (this->n_dofs(), false, type);
+          pr.second->init (this->n_dofs(), false, type);
         }
       else
         {
           libmesh_assert_equal_to(type, PARALLEL);
-          pos->second->init (this->n_dofs(), this->n_local_dofs(), false, type);
+          pr.second->init (this->n_dofs(), this->n_local_dofs(), false, type);
         }
     }
 }
@@ -323,30 +323,30 @@ void System::restrict_vectors ()
 {
 #ifdef LIBMESH_ENABLE_AMR
   // Restrict the _vectors on the coarsened cells
-  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
+  for (auto & pr : _vectors)
     {
-      NumericVector<Number> * v = pos->second;
+      NumericVector<Number> * v = pr.second;
 
-      if (_vector_projections[pos->first])
+      if (_vector_projections[pr.first])
         {
-          this->project_vector (*v, this->vector_is_adjoint(pos->first));
+          this->project_vector (*v, this->vector_is_adjoint(pr.first));
         }
       else
         {
-          ParallelType type = _vector_types[pos->first];
+          ParallelType type = _vector_types[pr.first];
 
           if (type == GHOSTED)
             {
 #ifdef LIBMESH_ENABLE_GHOSTED
-              pos->second->init (this->n_dofs(), this->n_local_dofs(),
-                                 _dof_map->get_send_list(), false,
-                                 GHOSTED);
+              pr.second->init (this->n_dofs(), this->n_local_dofs(),
+                               _dof_map->get_send_list(), false,
+                               GHOSTED);
 #else
               libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
 #endif
             }
           else
-            pos->second->init (this->n_dofs(), this->n_local_dofs(), false, type);
+            pr.second->init (this->n_dofs(), this->n_local_dofs(), false, type);
         }
     }
 
@@ -575,19 +575,18 @@ bool System::compare (const System & other_system,
   else
     {
       // compare other vectors
-      for (const_vectors_iterator pos = _vectors.begin();
-           pos != _vectors.end(); ++pos)
+      for (auto & pr : _vectors)
         {
           if (verbose)
             libMesh::out << "   comparing vector \""
-                         << pos->first << "\" ...";
+                         << pr.first << "\" ...";
 
           // assume they have the same name
           const NumericVector<Number> & other_system_vector =
-            other_system.get_vector(pos->first);
+            other_system.get_vector(pr.first);
 
-          ov_result.push_back(pos->second->compare (other_system_vector,
-                                                    threshold));
+          ov_result.push_back(pr.second->compare (other_system_vector,
+                                                  threshold));
 
           if (verbose)
             {
@@ -597,9 +596,7 @@ bool System::compare (const System & other_system,
                 libMesh::out << " first difference occurred at" << std::endl
                              << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
             }
-
         }
-
     } // finished comparing additional vectors
 
 
@@ -1483,18 +1480,17 @@ Real System::calculate_norm(const NumericVector<Number> & v,
       const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
 
       // Prepare finite elements for each dimension present in the mesh
-      for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
-           d_it != elem_dims.end(); ++d_it)
+      for (const auto & dim : elem_dims)
         {
-          if (skip_dimensions && skip_dimensions->find(*d_it) != skip_dimensions->end())
+          if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
             continue;
 
           // Construct quadrature and finite element objects
-          q_rules[*d_it] = fe_type.default_quadrature_rule (*d_it);
-          fe_ptrs[*d_it] = FEBase::build(*d_it, fe_type);
+          q_rules[dim] = fe_type.default_quadrature_rule (dim);
+          fe_ptrs[dim] = FEBase::build(dim, fe_type);
 
           // Attach quadrature rule to FE object
-          fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it].get());
+          fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
         }
 
       std::vector<dof_id_type> dof_indices;
diff --git a/src/systems/system_io.C b/src/systems/system_io.C
index 0710c06a70d..3e7417c20de 100644
--- a/src/systems/system_io.C
+++ b/src/systems/system_io.C
@@ -189,8 +189,8 @@ void System::read_header (Xdr & io,
             std::vector<subdomain_id_type> domain_array;
             if (this->processor_id() == 0)
               io.data (domain_array);
-            for (std::vector<subdomain_id_type>::iterator it = domain_array.begin(); it != domain_array.end(); ++it)
-              domains.insert(*it);
+            for (const auto & id : domain_array)
+              domains.insert(id);
           }
         this->comm().broadcast(domains);
 
@@ -572,25 +572,23 @@ void System::read_parallel_data (Xdr & io,
       if (this->variable(var).type().family != SCALAR)
         {
           // First read the node DOF values
-          for (std::vector<const DofObject *>::const_iterator
-                 it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
-            for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+          for (const auto & node : ordered_nodes)
+            for (unsigned int comp=0; comp<node->n_comp(sys_num, var); comp++)
               {
-                libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
                                              DofObject::invalid_id);
                 libmesh_assert_less (cnt, io_buffer.size());
-                this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
+                this->solution->set(node->dof_number(sys_num, var, comp), io_buffer[cnt++]);
               }
 
           // Then read the element DOF values
-          for (std::vector<const DofObject *>::const_iterator
-                 it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
-            for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+          for (const auto & elem : ordered_elements)
+            for (unsigned int comp=0; comp<elem->n_comp(sys_num, var); comp++)
               {
-                libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
                                              DofObject::invalid_id);
                 libmesh_assert_less (cnt, io_buffer.size());
-                this->solution->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
+                this->solution->set(elem->dof_number(sys_num, var, comp), io_buffer[cnt++]);
               }
         }
     }
@@ -661,25 +659,23 @@ void System::read_parallel_data (Xdr & io,
                   if (this->variable(var).type().family != SCALAR)
                     {
                       // First read the node DOF values
-                      for (std::vector<const DofObject *>::const_iterator
-                             it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
-                        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+                      for (const auto & node : ordered_nodes)
+                        for (unsigned int comp=0; comp<node->n_comp(sys_num, var); comp++)
                           {
-                            libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                            libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
                                                          DofObject::invalid_id);
                             libmesh_assert_less (cnt, io_buffer.size());
-                            pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
+                            pos->second->set(node->dof_number(sys_num, var, comp), io_buffer[cnt++]);
                           }
 
                       // Then read the element DOF values
-                      for (std::vector<const DofObject *>::const_iterator
-                             it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
-                        for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+                      for (const auto & elem : ordered_elements)
+                        for (unsigned int comp=0; comp<elem->n_comp(sys_num, var); comp++)
                           {
-                            libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                            libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
                                                          DofObject::invalid_id);
                             libmesh_assert_less (cnt, io_buffer.size());
-                            pos->second->set((*it)->dof_number(sys_num, var, comp), io_buffer[cnt++]);
+                            pos->second->set(elem->dof_number(sys_num, var, comp), io_buffer[cnt++]);
                           }
                     }
                 }
@@ -835,8 +831,6 @@ std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_obj
   // to write, and comps are all the components for said
   // vars on the object.
 
-  typedef std::vector<NumericVector<Number> *>::const_iterator vec_iterator_type;
-
   // variables to read.  Unless specified otherwise, defaults to _written_var_indices.
   std::vector<unsigned int> vars_to_read (_written_var_indices);
 
@@ -883,9 +877,8 @@ std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_obj
       xfer_ids_size[block] += 2; // for each object, we send its id, as well as the total number of components for all variables
 
       dof_id_type n_comp_tot=0;
-      for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
-           var_it!=vars_to_read.end(); ++var_it)
-        n_comp_tot += (*it)->n_comp(sys_num, *var_it); // for each variable, we will receive the nonzero components
+      for (const auto & var : vars_to_read)
+        n_comp_tot += (*it)->n_comp(sys_num, var); // for each variable, we will receive the nonzero components
 
       recv_vals_size[block] += n_comp_tot*num_vecs;
     }
@@ -926,9 +919,8 @@ std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_obj
 
               unsigned int n_comp_tot=0;
 
-              for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
-                   var_it!=vars_to_read.end(); ++var_it)
-                n_comp_tot += (*it)->n_comp(sys_num,*var_it);
+              for (const auto & var : vars_to_read)
+                n_comp_tot += (*it)->n_comp(sys_num, var);
 
               ids.push_back (n_comp_tot*num_vecs);
             }
@@ -1095,29 +1087,24 @@ std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_obj
             if (((*it)->id() >= first_object) && // object in [first_object,last_object)
                 ((*it)->id() <   last_object))
               // unpack & set the values
-              for (vec_iterator_type vec_it=vecs.begin(); vec_it!=vecs.end(); ++vec_it)
-                {
-                  NumericVector<Number> * vec(*vec_it);
-
-                  for (std::vector<unsigned int>::const_iterator var_it=vars_to_read.begin();
-                       var_it!=vars_to_read.end(); ++var_it)
-                    {
-                      const unsigned int n_comp = (*it)->n_comp(sys_num,*var_it);
+              for (auto & vec : vecs)
+                for (const auto & var : vars_to_read)
+                  {
+                    const unsigned int n_comp = (*it)->n_comp(sys_num, var);
 
-                      for (unsigned int comp=0; comp<n_comp; comp++, ++val_it)
-                        {
-                          const dof_id_type dof_index = (*it)->dof_number (sys_num, *var_it, comp);
-                          libmesh_assert (val_it != vals.end());
-                          if (vec)
-                            {
-                              libmesh_assert_greater_equal (dof_index, vec->first_local_index());
-                              libmesh_assert_less (dof_index, vec->last_local_index());
-                              //libMesh::out << "dof_index, *val_it = \t" << dof_index << ", " << *val_it << '\n';
-                              vec->set (dof_index, *val_it);
-                            }
-                        }
-                    }
-                }
+                    for (unsigned int comp=0; comp<n_comp; comp++, ++val_it)
+                      {
+                        const dof_id_type dof_index = (*it)->dof_number (sys_num, var, comp);
+                        libmesh_assert (val_it != vals.end());
+                        if (vec)
+                          {
+                            libmesh_assert_greater_equal (dof_index, vec->first_local_index());
+                            libmesh_assert_less (dof_index, vec->last_local_index());
+                            //libMesh::out << "dof_index, *val_it = \t" << dof_index << ", " << *val_it << '\n';
+                            vec->set (dof_index, *val_it);
+                          }
+                      }
+                  }
       }
 
       // processor 0 needs to make sure all replies have been handed off
@@ -1496,11 +1483,8 @@ void System::write_header (Xdr & io,
 
     if (write_additional_data)
       {
-        std::map<std::string, NumericVector<Number> *>::const_iterator
-          vec_pos = this->_vectors.begin();
         unsigned int cnt=0;
-
-        for (; vec_pos != this->_vectors.end(); ++vec_pos)
+        for (const auto & pr : _vectors)
           {
             // 9.)
             // write the name of the cnt-th additional vector
@@ -1508,7 +1492,7 @@ void System::write_header (Xdr & io,
             std::sprintf(buf, "%d", cnt++);
             comment += buf;
             comment += "th vector";
-            std::string vec_name = vec_pos->first;
+            std::string vec_name = pr.first;
 
             io.data (vec_name, comment.c_str());
           }
@@ -1586,26 +1570,23 @@ void System::write_parallel_data (Xdr & io,
     if (this->variable(var).type().family != SCALAR)
       {
         // First write the node DOF values
-        for (std::vector<const DofObject *>::const_iterator
-               it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
-          for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+        for (const auto & node : ordered_nodes)
+          for (unsigned int comp=0; comp<node->n_comp(sys_num, var); comp++)
             {
-              //libMesh::out << "(*it)->id()=" << (*it)->id() << std::endl;
-              libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+              libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
                                            DofObject::invalid_id);
 
-              io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
+              io_buffer.push_back((*this->solution)(node->dof_number(sys_num, var, comp)));
             }
 
         // Then write the element DOF values
-        for (std::vector<const DofObject *>::const_iterator
-               it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
-          for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+        for (const auto & elem : ordered_elements)
+          for (unsigned int comp=0; comp<elem->n_comp(sys_num, var); comp++)
             {
-              libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+              libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
                                            DofObject::invalid_id);
 
-              io_buffer.push_back((*this->solution)((*it)->dof_number(sys_num, var, comp)));
+              io_buffer.push_back((*this->solution)(elem->dof_number(sys_num, var, comp)));
             }
       }
 
@@ -1643,37 +1624,33 @@ void System::write_parallel_data (Xdr & io,
   // Only write additional vectors if wanted
   if (write_additional_data)
     {
-      std::map<std::string, NumericVector<Number> *>::const_iterator
-        pos = _vectors.begin();
-
-      for (; pos != this->_vectors.end(); ++pos)
+      for (auto & pr : _vectors)
         {
-          io_buffer.clear(); io_buffer.reserve( pos->second->local_size());
+          io_buffer.clear();
+          io_buffer.reserve(pr.second->local_size());
 
           // Loop over each non-SCALAR variable and each node, and write out the value.
           for (unsigned int var=0; var<nv; var++)
             if (this->variable(var).type().family != SCALAR)
               {
                 // First write the node DOF values
-                for (std::vector<const DofObject *>::const_iterator
-                       it = ordered_nodes.begin(); it != ordered_nodes.end(); ++it)
-                  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+                for (const auto & node : ordered_nodes)
+                  for (unsigned int comp=0; comp<node->n_comp(sys_num, var); comp++)
                     {
-                      libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                      libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
                                                    DofObject::invalid_id);
 
-                      io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
+                      io_buffer.push_back((*pr.second)(node->dof_number(sys_num, var, comp)));
                     }
 
                 // Then write the element DOF values
-                for (std::vector<const DofObject *>::const_iterator
-                       it = ordered_elements.begin(); it != ordered_elements.end(); ++it)
-                  for (unsigned int comp=0; comp<(*it)->n_comp(sys_num, var); comp++)
+                for (const auto & elem : ordered_elements)
+                  for (unsigned int comp=0; comp<elem->n_comp(sys_num, var); comp++)
                     {
-                      libmesh_assert_not_equal_to ((*it)->dof_number(sys_num, var, comp),
+                      libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
                                                    DofObject::invalid_id);
 
-                      io_buffer.push_back((*pos->second)((*it)->dof_number(sys_num, var, comp)));
+                      io_buffer.push_back((*pr.second)(elem->dof_number(sys_num, var, comp)));
                     }
               }
 
@@ -1688,7 +1665,7 @@ void System::write_parallel_data (Xdr & io,
                     dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
 
                     for (std::size_t i=0; i<SCALAR_dofs.size(); i++)
-                      io_buffer.push_back((*pos->second)(SCALAR_dofs[i]));
+                      io_buffer.push_back((*pr.second)(SCALAR_dofs[i]));
                   }
               }
 
@@ -1702,7 +1679,7 @@ void System::write_parallel_data (Xdr & io,
             comment = "# System \"";
             comment += this->name();
             comment += "\" Additional Vector \"";
-            comment += pos->first;
+            comment += pr.first;
             comment += "\"";
           }
 
@@ -1858,8 +1835,6 @@ std::size_t System::write_serialized_blocked_dof_objects (const std::vector<cons
   // to write, and comps are all the components for said
   // vars on the object.
 
-  typedef std::vector<const NumericVector<Number> *>::const_iterator vec_iterator_type;
-
   // We will write all variables unless requested otherwise.
   std::vector<unsigned int> vars_to_write(1, var_to_write);
 
@@ -1910,9 +1885,8 @@ std::size_t System::write_serialized_blocked_dof_objects (const std::vector<cons
 
       unsigned int n_comp_tot=0;
 
-      for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
-           var_it!=vars_to_write.end(); ++var_it)
-        n_comp_tot += (*it)->n_comp(sys_num, *var_it); // for each variable, we will store the nonzero components
+      for (const auto & var : vars_to_write)
+        n_comp_tot += (*it)->n_comp(sys_num, var); // for each variable, we will store the nonzero components
 
       send_vals_size[block] += n_comp_tot*num_vecs;
     }
@@ -1951,31 +1925,25 @@ std::size_t System::write_serialized_blocked_dof_objects (const std::vector<cons
               {
                 unsigned int n_comp_tot=0;
 
-                for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
-                     var_it!=vars_to_write.end(); ++var_it)
-                  n_comp_tot += (*it)->n_comp(sys_num,*var_it);
+                for (const auto & var : vars_to_write)
+                  n_comp_tot += (*it)->n_comp(sys_num, var);
 
                 ids.push_back (n_comp_tot*num_vecs); // even if 0 - processor 0 has no way of knowing otherwise...
               }
 
               // pack the values to send
-              for (vec_iterator_type vec_it=vecs.begin(); vec_it!=vecs.end(); ++vec_it)
-                {
-                  const NumericVector<Number> & vec(**vec_it);
-
-                  for (std::vector<unsigned int>::const_iterator var_it=vars_to_write.begin();
-                       var_it!=vars_to_write.end(); ++var_it)
-                    {
-                      const unsigned int n_comp = (*it)->n_comp(sys_num,*var_it);
+              for (const auto & vec : vecs)
+                for (const auto & var : vars_to_write)
+                  {
+                    const unsigned int n_comp = (*it)->n_comp(sys_num, var);
 
-                      for (unsigned int comp=0; comp<n_comp; comp++)
-                        {
-                          libmesh_assert_greater_equal ((*it)->dof_number(sys_num, *var_it, comp), vec.first_local_index());
-                          libmesh_assert_less ((*it)->dof_number(sys_num, *var_it, comp), vec.last_local_index());
-                          vals.push_back(vec((*it)->dof_number(sys_num, *var_it, comp)));
-                        }
-                    }
-                }
+                    for (unsigned int comp=0; comp<n_comp; comp++)
+                      {
+                        libmesh_assert_greater_equal ((*it)->dof_number(sys_num, var, comp), vec->first_local_index());
+                        libmesh_assert_less ((*it)->dof_number(sys_num, var, comp), vec->last_local_index());
+                        vals.push_back((*vec)((*it)->dof_number(sys_num, var, comp)));
+                      }
+                  }
             }
 
 #ifdef LIBMESH_HAVE_MPI
diff --git a/src/systems/system_projection.C b/src/systems/system_projection.C
index 748eff642c0..a0b011f66f5 100644
--- a/src/systems/system_projection.C
+++ b/src/systems/system_projection.C
@@ -1086,9 +1086,8 @@ void BuildProjectionList::operator()(const ConstElemRange & range)
   std::vector<dof_id_type> di;
 
   // Iterate over the elements in the range
-  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
+  for (const auto & elem : range)
     {
-      const Elem * elem = *elem_it;
       // If this element doesn't have an old_dof_object with dofs for the
       // current system, then it must be newly added, so the user
       // is responsible for setting the new dofs.
@@ -1266,10 +1265,8 @@ void BoundaryProjectSolution::operator()(const ConstElemRange & range) const
       std::vector<boundary_id_type> bc_ids;
 
       // Iterate over all the elements in the range
-      for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
+      for (const auto & elem : range)
         {
-          const Elem * elem = *elem_it;
-
           // Per-subdomain variables don't need to be projected on
           // elements where they're not active
           if (!variable.active_on_subdomain(elem->subdomain_id()))
@@ -1289,9 +1286,8 @@ void BoundaryProjectSolution::operator()(const ConstElemRange & range) const
               // First see if this side has been requested
               boundary_info.boundary_ids (elem, s, bc_ids);
               bool do_this_side = false;
-              for (std::vector<boundary_id_type>::iterator i=bc_ids.begin();
-                   i!=bc_ids.end(); ++i)
-                if (b.count(*i))
+              for (const auto & bc_id : bc_ids)
+                if (b.count(bc_id))
                   {
                     do_this_side = true;
                     break;
diff --git a/src/utils/point_locator_base.C b/src/utils/point_locator_base.C
index a707aa9e8d9..e942991d20d 100644
--- a/src/utils/point_locator_base.C
+++ b/src/utils/point_locator_base.C
@@ -105,11 +105,8 @@ locate_node(const Point & p,
   std::set<const Elem *> candidate_elements;
   this->operator()(p, candidate_elements, allowed_subdomains);
 
-  for (std::set<const Elem *>::const_iterator
-         it = candidate_elements.begin();
-       it != candidate_elements.end(); ++it)
+  for (const auto & elem : candidate_elements)
     {
-      const Elem * elem = *it;
       const int elem_n_nodes = elem->n_nodes();
       const Real hmax = elem->hmax();
       const Real dist_tol_sq = (tol * hmax) * (tol * hmax);
diff --git a/src/utils/tree_node.C b/src/utils/tree_node.C
index b2dd1eff39f..85300498f60 100644
--- a/src/utils/tree_node.C
+++ b/src/utils/tree_node.C
@@ -370,9 +370,8 @@ void TreeNode<N>::print_elements(std::ostream & out_stream) const
     {
       out_stream << "TreeNode Level: " << this->level() << std::endl;
 
-      for (std::vector<const Elem *>::const_iterator pos=elements.begin();
-           pos != elements.end(); ++pos)
-        out_stream << " " << *pos;
+      for (const auto & elem : elements)
+        out_stream << " " << elem;
 
       out_stream << std::endl << std::endl;
     }
@@ -419,16 +418,15 @@ void TreeNode<N>::transform_nodes_to_elements (std::vector<std::vector<const Ele
       // than the set.
       elements.reserve(elements_set.size());
 
-      for (std::set<const Elem *>::iterator pos=elements_set.begin();
-           pos != elements_set.end(); ++pos)
+      for (const auto & elem : elements_set)
         {
-          elements.push_back(*pos);
+          elements.push_back(elem);
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
           // flag indicating this node contains
           // infinite elements
-          if ((*pos)->infinite())
+          if (elem->infinite())
             this->contains_ifems = true;
 
 #endif
@@ -475,11 +473,10 @@ TreeNode<N>::find_element (const Point & p,
       // or if the node contains infinite elements
       if (this->bounds_point(p, relative_tol) || this->contains_ifems)
         // Search the active elements in the active TreeNode.
-        for (std::vector<const Elem *>::const_iterator pos=elements.begin();
-             pos != elements.end(); ++pos)
-          if (!allowed_subdomains || allowed_subdomains->count((*pos)->subdomain_id()))
-            if ((*pos)->active() && (*pos)->contains_point(p, relative_tol))
-              return *pos;
+        for (const auto & elem : elements)
+          if (!allowed_subdomains || allowed_subdomains->count(elem->subdomain_id()))
+            if (elem->active() && elem->contains_point(p, relative_tol))
+              return elem;
 
       // The point was not found in any element
       return libmesh_nullptr;