diff --git a/src/ASM/BlockElmMats.C b/src/ASM/BlockElmMats.C
index 9dd8eb8e..6e3fa26d 100644
--- a/src/ASM/BlockElmMats.C
+++ b/src/ASM/BlockElmMats.C
@@ -73,20 +73,14 @@ bool BlockElmMats::redimOffDiag (size_t blkIndex, char symmetric)
 }
 
 
-bool BlockElmMats::redimNewtonMat ()
+bool BlockElmMats::finalize ()
 {
   // Calculate the total number of element dofs
-  size_t neldof = 0;
+  neldof = 0;
   size_t nDiagB = blockInfo.size();
   for (size_t i = 0; i < nDiagB; i++)
     neldof += blockInfo[i].ncmp*basisInfo[blockInfo[i].basis-1].nen;
 
-  // Set the element Newton matrix dimension
-  if (!A.empty() && A.front().empty())
-    A.front().resize(neldof,neldof);
-  if (!b.empty() && b.front().empty())
-    b.front().resize(neldof);
-
   // Calculate the offset of each block sub-matrix
   size_t idof = 1;
   for (size_t j = 0; j < basisInfo.size(); j++)
@@ -103,13 +97,18 @@ bool BlockElmMats::redimNewtonMat ()
     idof += ndof;
   }
 
-  return A.empty() ? true : idof-1 == A.front().rows();
+  return idof-1 == neldof;
 }
 
 
 const Matrix& BlockElmMats::getNewtonMatrix () const
 {
+  if (!A.front().empty())
+    return A.front();
+
   Matrix& N = const_cast<Matrix&>(A.front());
+  // Set the element Newton matrix dimension
+  N.resize(neldof,neldof);
 
   size_t ib, jb, kb;
   size_t nDiagB = blockInfo.size();
@@ -174,7 +173,11 @@ const Matrix& BlockElmMats::getNewtonMatrix () const
 
 const Vector& BlockElmMats::getRHSVector () const
 {
+  if (!b.front().empty())
+    return b.front();
+
   Vector& R = const_cast<Vector&>(b.front());
+  R.resize(neldof);
 
   for (size_t ib = 0; ib < blockInfo.size(); ib++)
   {
diff --git a/src/ASM/BlockElmMats.h b/src/ASM/BlockElmMats.h
index 2298f516..0e7dcb2f 100644
--- a/src/ASM/BlockElmMats.h
+++ b/src/ASM/BlockElmMats.h
@@ -32,7 +32,8 @@ public:
   //! \brief The constructor initializes the block information arrays.
   //! \param[in] nBlk Number of matrix blocks (in each direction, row & column)
   //! \param[in] nBas Number of bases (> 1 for mixed problems)
-  BlockElmMats(size_t nBlk, size_t nBas) : blockInfo(nBlk), basisInfo(nBas) {}
+  BlockElmMats(size_t nBlk, size_t nBas) :
+    neldof(0), blockInfo(nBlk), basisInfo(nBas) {}
   //! \brief Empty destructor.
   virtual ~BlockElmMats() {}
 
@@ -52,10 +53,10 @@ public:
   //! \details This method must not be invoked until the \a redim method has
   //! been invoked for all the diagonal blocks.
   bool redimOffDiag(size_t blkIndex, char symmetric = 1);
-  //! \brief Sets the dimension of the Newton matrix.
+  //! \brief Calculates the dimension of the Newton matrix.
   //! \details This method must not be invoked until the \a redim method has
   //! been invoked for all the diagonal blocks.
-  bool redimNewtonMat();
+  bool finalize();
 
   //! \brief Returns the element-level Newton matrix.
   virtual const Matrix& getNewtonMatrix() const;
@@ -84,6 +85,7 @@ private:
     Basis() : nen(0), ncmp(0) {}
   };
 
+  size_t neldof; //!< Size of newton matrix
   std::vector<Block> blockInfo; //!< Block information
   std::vector<Basis> basisInfo; //!< Basis information
   std::vector<char>  symmFlag;  //!< Symmetry flag for off-diagonal blocks
diff --git a/src/ASM/Test/TestBlockElmMats.C b/src/ASM/Test/TestBlockElmMats.C
index e79a8d0a..18676cf9 100644
--- a/src/ASM/Test/TestBlockElmMats.C
+++ b/src/ASM/Test/TestBlockElmMats.C
@@ -22,7 +22,7 @@ TEST(TestBlockElmMats, 1Basis2BlocksDiag)
   mats.resize(3, 3);
   ASSERT_TRUE(mats.redim(1, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -50,7 +50,7 @@ TEST(TestBlockElmMats, 1Basis2BlocksSymmetric)
   ASSERT_TRUE(mats.redim(1, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1));
   ASSERT_TRUE(mats.redimOffDiag(3, 1));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -82,7 +82,7 @@ TEST(TestBlockElmMats, 1Basis2BlocksSkewSymmetric)
   ASSERT_TRUE(mats.redim(1, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1));
   ASSERT_TRUE(mats.redimOffDiag(3, -1));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -115,7 +115,7 @@ TEST(TestBlockElmMats, 1Basis2BlocksFull)
   ASSERT_TRUE(mats.redim(2, 2, 1));
   ASSERT_TRUE(mats.redimOffDiag(3, 0));
   ASSERT_TRUE(mats.redimOffDiag(4, 0));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -146,7 +146,7 @@ TEST(TestBlockElmMats, 2Basis2BlocksDiag)
   mats.resize(3, 3);
   ASSERT_TRUE(mats.redim(1, 2, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1, 2));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -173,7 +173,7 @@ TEST(TestBlockElmMats, 2Basis2BlocksSymmetric)
   ASSERT_TRUE(mats.redim(1, 2, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1, 2));
   ASSERT_TRUE(mats.redimOffDiag(3, 1));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -201,7 +201,7 @@ TEST(TestBlockElmMats, 2Basis2BlocksSkewSymmetric)
   ASSERT_TRUE(mats.redim(1, 2, 2, 1));
   ASSERT_TRUE(mats.redim(2, 2, 1, 2));
   ASSERT_TRUE(mats.redimOffDiag(3, -1));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);
@@ -230,7 +230,7 @@ TEST(TestBlockElmMats, 2Basis2BlocksFull)
   ASSERT_TRUE(mats.redim(2, 2, 1, 2));
   ASSERT_TRUE(mats.redimOffDiag(3, 0));
   ASSERT_TRUE(mats.redimOffDiag(4, 0));
-  mats.redimNewtonMat();
+  mats.finalize();
 
   mats.A[1].fill(1);
   mats.A[2].fill(2);