diff --git a/include/numerics/petsc_vector.h b/include/numerics/petsc_vector.h
index bd3ea1f5a6f..f8e53075edb 100644
--- a/include/numerics/petsc_vector.h
+++ b/include/numerics/petsc_vector.h
@@ -293,6 +293,30 @@ class PetscVector libmesh_final : public NumericVector<T>
   virtual void get(const std::vector<numeric_index_type> & index,
                    T * values) const libmesh_override;
 
+  /**
+   * Get read/write access to the raw PETSc Vector data array
+   *
+   * Note: This is an extremely advanced interface.
+   * In general you should avoid using it!
+   */
+  PetscScalar * get_array();
+
+  /**
+   * Get read only access to the raw PETSc Vector data array
+   *
+   * Note: This is an extremely advanced interface.
+   * In general you should avoid using it!
+   */
+  const PetscScalar * get_array_read() const;
+
+  /**
+   * Restore the data array.
+   *
+   * MUST be called after get_array() or get_array_read() and before using
+   * any other interface functions on PetscVector.
+   */
+  void restore_array();
+
   /**
    * Addition operator.
    * Fast equivalent to \p U.add(1, V).
@@ -540,14 +564,10 @@ class PetscVector libmesh_final : public NumericVector<T>
    */
   mutable numeric_index_type _last;
 
-#ifndef NDEBUG
   /**
-   * Size of the local form, for being used in assertations.  The
-   * contents of this field are only valid if the vector is ghosted
-   * and \p _array_is_present is \p true.
+   * Size of the local values from _get_array()
    */
   mutable numeric_index_type _local_size;
-#endif
 
   /**
    * Petsc vector datatype to hold the local form of a ghosted vector.
@@ -564,7 +584,17 @@ class PetscVector libmesh_final : public NumericVector<T>
    * const member functions, so _values also needs to be mutable
    * (otherwise it is a "const PetscScalar * const" in that context).
    */
-  mutable const PetscScalar * _values;
+  mutable const PetscScalar * _read_only_values;
+
+  /**
+   * Pointer to the actual Petsc array of the values of the vector.
+   * This pointer is only valid if \p _array_is_present is \p true.
+   * We're using Petsc's VecGetArrayRead() function, which requires a
+   * constant PetscScalar *, but _get_array and _restore_array are
+   * const member functions, so _values also needs to be mutable
+   * (otherwise it is a "const PetscScalar * const" in that context).
+   */
+  mutable PetscScalar * _values;
 
   /**
    * Mutex for _get_array and _restore_array.  This is part of the
@@ -580,8 +610,10 @@ class PetscVector libmesh_final : public NumericVector<T>
   /**
    * Queries the array (and the local form if the vector is ghosted)
    * from Petsc.
+   *
+   * @param read_only Whether or not a read only copy of the raw data
    */
-  void _get_array() const;
+  void _get_array(bool read_only) const;
 
   /**
    * Restores the array (and the local form if the vector is ghosted)
@@ -605,6 +637,17 @@ class PetscVector libmesh_final : public NumericVector<T>
    * for the constructor which takes a PETSc Vec object.
    */
   bool _destroy_vec_on_exit;
+
+  /**
+   * Whether or not the data array has been manually retrieved using
+   * get_array() or get_array_read()
+   */
+  mutable bool _values_manually_retrieved;
+
+  /**
+   * Whether or not the data array is for read only access
+   */
+  mutable bool _values_read_only;
 };
 
 
@@ -622,7 +665,9 @@ PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in, const Paral
   _local_form(libmesh_nullptr),
   _values(libmesh_nullptr),
   _global_to_local_map(),
-  _destroy_vec_on_exit(true)
+  _destroy_vec_on_exit(true),
+  _values_manually_retrieved(false),
+  _values_read_only(false)
 {
   this->_type = ptype;
 }
@@ -639,7 +684,9 @@ PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
   _local_form(libmesh_nullptr),
   _values(libmesh_nullptr),
   _global_to_local_map(),
-  _destroy_vec_on_exit(true)
+  _destroy_vec_on_exit(true),
+  _values_manually_retrieved(false),
+  _values_read_only(false)
 {
   this->init(n, n, false, ptype);
 }
@@ -657,7 +704,9 @@ PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
   _local_form(libmesh_nullptr),
   _values(libmesh_nullptr),
   _global_to_local_map(),
-  _destroy_vec_on_exit(true)
+  _destroy_vec_on_exit(true),
+  _values_manually_retrieved(false),
+  _values_read_only(false)
 {
   this->init(n, n_local, false, ptype);
 }
@@ -676,7 +725,9 @@ PetscVector<T>::PetscVector (const Parallel::Communicator & comm_in,
   _local_form(libmesh_nullptr),
   _values(libmesh_nullptr),
   _global_to_local_map(),
-  _destroy_vec_on_exit(true)
+  _destroy_vec_on_exit(true),
+  _values_manually_retrieved(false),
+  _values_read_only(false)
 {
   this->init(n, n_local, ghost, false, ptype);
 }
@@ -694,7 +745,9 @@ PetscVector<T>::PetscVector (Vec v,
   _local_form(libmesh_nullptr),
   _values(libmesh_nullptr),
   _global_to_local_map(),
-  _destroy_vec_on_exit(false)
+  _destroy_vec_on_exit(false),
+  _values_manually_retrieved(false),
+  _values_read_only(false)
 {
   this->_vec = v;
   this->_is_closed = true;
@@ -1219,7 +1272,7 @@ template <typename T>
 inline
 T PetscVector<T>::operator() (const numeric_index_type i) const
 {
-  this->_get_array();
+  this->_get_array(true);
 
   const numeric_index_type local_index = this->map_global_to_local_index(i);
 
@@ -1230,7 +1283,7 @@ T PetscVector<T>::operator() (const numeric_index_type i) const
     }
 #endif
 
-  return static_cast<T>(_values[local_index]);
+  return static_cast<T>(_read_only_values[local_index]);
 }
 
 
@@ -1240,7 +1293,7 @@ inline
 void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
                          T * values) const
 {
-  this->_get_array();
+  this->_get_array(true);
 
   const std::size_t num = index.size();
 
@@ -1253,11 +1306,41 @@ void PetscVector<T>::get(const std::vector<numeric_index_type> & index,
           libmesh_assert_less (local_index, _local_size);
         }
 #endif
-      values[i] = static_cast<T>(_values[local_index]);
+      values[i] = static_cast<T>(_read_only_values[local_index]);
     }
 }
 
 
+template <typename T>
+inline
+PetscScalar * PetscVector<T>::get_array()
+{
+  _values_manually_retrieved = true;
+  _get_array(false);
+
+  return _values;
+}
+
+
+template <typename T>
+inline
+const PetscScalar * PetscVector<T>::get_array_read() const
+{
+  _values_manually_retrieved = true;
+  _get_array(true);
+
+  return _read_only_values;
+}
+
+template <typename T>
+inline
+void PetscVector<T>::restore_array()
+{
+  // Note _values_manually_retrieved needs to be set to false
+  // BEFORE calling _restore_array()!
+  _values_manually_retrieved = false;
+  _restore_array();
+}
 
 template <typename T>
 inline
diff --git a/src/numerics/petsc_vector.C b/src/numerics/petsc_vector.C
index ac86a963398..7d7ba21128d 100644
--- a/src/numerics/petsc_vector.C
+++ b/src/numerics/petsc_vector.C
@@ -299,64 +299,10 @@ void PetscVector<T>::add_vector_conjugate_transpose (const NumericVector<T> & V_
 template <typename T>
 void PetscVector<T>::add (const T v_in)
 {
-  this->_restore_array();
-
-  PetscErrorCode ierr=0;
-  PetscScalar * values;
-  const PetscScalar v = static_cast<PetscScalar>(v_in);
-
-  if(this->type() != GHOSTED)
-    {
-      const PetscInt n   = static_cast<PetscInt>(this->local_size());
-      const PetscInt fli = static_cast<PetscInt>(this->first_local_index());
-
-      for (PetscInt i=0; i<n; i++)
-        {
-          ierr = VecGetArray (_vec, &values);
-          LIBMESH_CHKERR(ierr);
-
-          PetscInt ig = fli + i;
-
-          PetscScalar value = (values[i] + v);
-
-          ierr = VecRestoreArray (_vec, &values);
-          LIBMESH_CHKERR(ierr);
-
-          ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
-          LIBMESH_CHKERR(ierr);
-        }
-    }
-  else
-    {
-      /* Vectors that include ghost values require a special
-         handling.  */
-      Vec loc_vec;
-      ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-      LIBMESH_CHKERR(ierr);
-
-      PetscInt n=0;
-      ierr = VecGetSize(loc_vec, &n);
-      LIBMESH_CHKERR(ierr);
-
-      for (PetscInt i=0; i<n; i++)
-        {
-          ierr = VecGetArray (loc_vec, &values);
-          LIBMESH_CHKERR(ierr);
-
-          PetscScalar value = (values[i] + v);
-
-          ierr = VecRestoreArray (loc_vec, &values);
-          LIBMESH_CHKERR(ierr);
-
-          ierr = VecSetValues (loc_vec, 1, &i, &value, INSERT_VALUES);
-          LIBMESH_CHKERR(ierr);
-        }
-
-      ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-      LIBMESH_CHKERR(ierr);
-    }
+  this->_get_array(false);
 
-  this->_is_closed = false;
+  for (numeric_index_type i=0; i<_local_size; i++)
+    _values[i] += v_in;
 }
 
 
@@ -655,12 +601,7 @@ template <typename T>
 NumericVector<T> &
 PetscVector<T>::operator = (const std::vector<T> & v)
 {
-  this->_restore_array();
-
-  const numeric_index_type nl   = this->local_size();
-  const numeric_index_type ioff = this->first_local_index();
-  PetscErrorCode ierr=0;
-  PetscScalar * values;
+  this->_get_array(false);
 
   /**
    * Case 1:  The vector is the same size of
@@ -668,14 +609,8 @@ PetscVector<T>::operator = (const std::vector<T> & v)
    */
   if (this->size() == v.size())
     {
-      ierr = VecGetArray (_vec, &values);
-      LIBMESH_CHKERR(ierr);
-
-      for (numeric_index_type i=0; i<nl; i++)
-        values[i] =  static_cast<PetscScalar>(v[i+ioff]);
-
-      ierr = VecRestoreArray (_vec, &values);
-      LIBMESH_CHKERR(ierr);
+      for (numeric_index_type i=0; i<_local_size; i++)
+        _values[i] =  static_cast<PetscScalar>(v[_first + i]);
     }
 
   /**
@@ -684,16 +619,8 @@ PetscVector<T>::operator = (const std::vector<T> & v)
    */
   else
     {
-      libmesh_assert_equal_to (this->local_size(), v.size());
-
-      ierr = VecGetArray (_vec, &values);
-      LIBMESH_CHKERR(ierr);
-
-      for (numeric_index_type i=0; i<nl; i++)
-        values[i] = static_cast<PetscScalar>(v[i]);
-
-      ierr = VecRestoreArray (_vec, &values);
-      LIBMESH_CHKERR(ierr);
+      for (numeric_index_type i=0; i<_local_size; i++)
+        _values[i] = static_cast<PetscScalar>(v[i]);
     }
 
   // Make sure ghost dofs are up to date
@@ -1395,7 +1322,7 @@ void PetscVector<T>::create_subvector(NumericVector<T> & subvector,
 
 
 template <typename T>
-void PetscVector<T>::_get_array() const
+void PetscVector<T>::_get_array(bool read_only) const
 {
   libmesh_assert (this->initialized());
 
@@ -1405,6 +1332,18 @@ void PetscVector<T>::_get_array() const
   Threads::spin_mutex::scoped_lock lock(_petsc_vector_mutex);
 #endif
 
+  // If the values have already been retrieved and we're currently
+  // trying to get a non-read only view (ie read/write) and the
+  // values are currently read only... then we need to restore
+  // the array first... and then retrieve it again.
+  if (_array_is_present && !read_only && _values_read_only)
+    _restore_array();
+
+  // If we already have a read/write array - and we're trying
+  // to get a read only array - let's just use the read write
+  if (_array_is_present && read_only && !_values_read_only)
+    _read_only_values = _values;
+
   if (!_array_is_present)
     {
 #ifdef LIBMESH_HAVE_CXX11_THREAD
@@ -1420,10 +1359,21 @@ void PetscVector<T>::_get_array() const
               // have an older PETSc than that, we'll do an ugly
               // const_cast and call VecGetArray() instead.
               ierr = VecGetArray(_vec, const_cast<PetscScalar **>(&_values));
+              _values_read_only = false;
 #else
-              ierr = VecGetArrayRead(_vec, &_values);
+              if (read_only)
+              {
+                ierr = VecGetArrayRead(_vec, &_read_only_values);
+                _values_read_only = true;
+              }
+              else
+              {
+                ierr = VecGetArray(_vec, &_values);
+                _values_read_only = false;
+              }
 #endif
               LIBMESH_CHKERR(ierr);
+              _local_size = this->local_size();
             }
           else
             {
@@ -1435,16 +1385,25 @@ void PetscVector<T>::_get_array() const
               // have an older PETSc than that, we'll do an ugly
               // const_cast and call VecGetArray() instead.
               ierr = VecGetArray(_local_form, const_cast<PetscScalar **>(&_values));
+              _values_read_only = false;
 #else
-              ierr = VecGetArrayRead(_local_form, &_values);
+              if (read_only)
+              {
+                ierr = VecGetArrayRead(_local_form, &_read_only_values);
+                _values_read_only = true;
+              }
+              else
+              {
+                ierr = VecGetArray(_local_form, &_values);
+                _values_read_only = false;
+              }
 #endif
               LIBMESH_CHKERR(ierr);
-#ifndef NDEBUG
+
               PetscInt my_local_size = 0;
               ierr = VecGetLocalSize(_local_form, &my_local_size);
               LIBMESH_CHKERR(ierr);
               _local_size = static_cast<numeric_index_type>(my_local_size);
-#endif
             }
 
           { // cache ownership range
@@ -1469,6 +1428,9 @@ void PetscVector<T>::_get_array() const
 template <typename T>
 void PetscVector<T>::_restore_array() const
 {
+  if (_values_manually_retrieved)
+    libmesh_error_msg("PetscVector values were manually retrieved but are being automatically restored!");
+
 #ifdef LIBMESH_HAVE_CXX11_THREAD
   std::atomic_thread_fence(std::memory_order_acquire);
 #else
@@ -1493,7 +1455,10 @@ void PetscVector<T>::_restore_array() const
               // const_cast and call VecRestoreArray() instead.
               ierr = VecRestoreArray (_vec, const_cast<PetscScalar **>(&_values));
 #else
-              ierr = VecRestoreArrayRead (_vec, &_values);
+              if (_values_read_only)
+                ierr = VecRestoreArrayRead (_vec, &_read_only_values);
+              else
+                ierr = VecRestoreArray (_vec, &_values);
 #endif
 
               LIBMESH_CHKERR(ierr);
@@ -1507,16 +1472,17 @@ void PetscVector<T>::_restore_array() const
               // const_cast and call VecRestoreArray() instead.
               ierr = VecRestoreArray (_local_form, const_cast<PetscScalar **>(&_values));
 #else
-              ierr = VecRestoreArrayRead (_local_form, &_values);
+              if (_values_read_only)
+                ierr = VecRestoreArrayRead (_local_form, &_read_only_values);
+              else
+                ierr = VecRestoreArray (_local_form, &_values);
 #endif
               LIBMESH_CHKERR(ierr);
               _values = libmesh_nullptr;
               ierr = VecGhostRestoreLocalForm (_vec,&_local_form);
               LIBMESH_CHKERR(ierr);
               _local_form = libmesh_nullptr;
-#ifndef NDEBUG
               _local_size = 0;
-#endif
             }
 #ifdef LIBMESH_HAVE_CXX11_THREAD
           std::atomic_thread_fence(std::memory_order_release);
diff --git a/tests/numerics/petsc_vector_test.C b/tests/numerics/petsc_vector_test.C
index 6ac097aa619..fe9c0ab582a 100644
--- a/tests/numerics/petsc_vector_test.C
+++ b/tests/numerics/petsc_vector_test.C
@@ -27,7 +27,52 @@ public:
 
   NUMERICVECTORTEST
 
+  CPPUNIT_TEST( testGetArray );
+
   CPPUNIT_TEST_SUITE_END();
+
+  void testGetArray()
+  {
+    unsigned int block_size  = 2;
+
+    // a different size on each processor.
+    unsigned int local_size  = block_size;
+    unsigned int global_size = 0;
+
+    for (libMesh::processor_id_type p=0; p<my_comm->size(); p++)
+      global_size += (block_size + static_cast<unsigned int>(p));
+
+    PetscVector<Number> v(*my_comm, global_size, local_size);
+
+    PetscScalar * values = v.get_array();
+
+    for (unsigned int i=0; i<local_size; i++)
+      values[i] = i;
+
+    v.restore_array();
+
+    v.close();
+
+    // Check the values through the interface
+    for (unsigned int i=0; i<local_size; i++)
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(v(my_comm->rank()*2 + i), i, TOLERANCE*TOLERANCE);
+
+    // Check that we can see the same thing with get_array_read
+    const PetscScalar * read_only_values = v.get_array_read();
+
+    for (unsigned int i=0; i<local_size; i++)
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(read_only_values[i], i, TOLERANCE*TOLERANCE);
+
+    v.restore_array();
+
+    // Test getting a read only array after getting a writeable array
+    values = v.get_array();
+    read_only_values = v.get_array_read();
+    CPPUNIT_ASSERT_EQUAL((intptr_t)read_only_values, (intptr_t)values);
+
+    v.restore_array();
+  }
+
 };
 
 CPPUNIT_TEST_SUITE_REGISTRATION( PetscVectorTest );