diff --git a/Apps/Common/CompatibleOperators.C b/Apps/Common/CompatibleOperators.C
index 02c8d477..c0f160ea 100644
--- a/Apps/Common/CompatibleOperators.C
+++ b/Apps/Common/CompatibleOperators.C
@@ -22,10 +22,9 @@ void CompatibleOperators::Weak::Advection(std::vector<Matrix>& EM,
 {
   size_t nsd = fe.grad(1).cols();
   for (size_t n = 1; n <= nsd; ++n)
-    for (size_t i = 1; i <= fe.basis(n).size(); ++i)
-      for (size_t j = 1; j <= fe.basis(n).size(); ++j)
-        for (size_t k = 1; k <= nsd; ++k)
-          EM[n](i,j) += scale*AC[k]*fe.grad(n)(j,k)*fe.basis(n)(i)*fe.detJxW;
+    for (size_t k = 1; k <= nsd; ++k)
+        EM[n].outer_product(fe.basis(n), fe.grad(n).getColumn(k),
+                            true,scale*AC[k]*fe.detJxW);
 }
 
 
@@ -59,9 +58,8 @@ void CompatibleOperators::Weak::Gradient(std::vector<Matrix>& EM,
 {
   size_t nsd = fe.grad(1).cols();
   for (size_t n = 1; n <= nsd; ++n)
-    for (size_t i=1; i <= fe.basis(n).size(); ++i)
-      for (size_t j=1; j <= fe.basis(nsd+1).size(); ++j)
-        EM[10+5*(n-1)](i,j) += -scale*fe.grad(n)(i,n)*fe.basis(nsd+1)(j)*fe.detJxW;
+    EM[10+5*(n-1)].outer_product(fe.grad(n).getColumn(n), fe.basis(nsd+1),
+                                 true, -scale*fe.detJxW);
 }
 
 
@@ -69,8 +67,7 @@ void CompatibleOperators::Weak::Gradient(Vectors& EV, const FiniteElement& fe,
                                          double scale)
 {
   for (size_t n = 1; n <= fe.grad(n).cols(); ++n)
-    for (size_t i = 1; i <= fe.basis(n).size(); ++i)
-      EV[n](i) += scale*fe.grad(n)(i,n)*fe.detJxW;
+    EV[n].add(fe.grad(n).getColumn(n), scale*fe.detJxW);
 }
 
 
@@ -83,12 +80,11 @@ void CompatibleOperators::Weak::Laplacian(std::vector<Matrix>& EM,
     EqualOrderOperators::Weak::Laplacian(EM[n], fe, scale, false, n);
 
   for (size_t m = 1; m <= nsd && stress; m++)
-    for (size_t i = 1; i <= fe.basis(m).size(); ++i)
-      for (size_t n = m; n <= nsd; n++)
-        for (size_t j = 1; j <= fe.basis(n).size(); ++j) {
-          int idx = m == n ? m : (m == 1 ? 6+n-m : 13);
-          EM[idx](i,j) += scale* fe.grad(n)(j,m) * fe.grad(m)(i,n) * fe.detJxW;
-        }
+    for (size_t n = m; n <= nsd; n++) {
+      int idx = m == n ? m : (m == 1 ? 6+n-m : 13);
+      EM[idx].outer_product(fe.grad(m).getColumn(n),
+                            fe.grad(n).getColumn(m), true, scale*fe.detJxW);
+    }
 }
 
 
@@ -127,8 +123,7 @@ void CompatibleOperators::Residual::Convection(Vectors& EV, const FiniteElement&
     double conv = 0.0;
     for (size_t k = 1; k<= nsd; ++k)
       conv += U[k-1]*dUdX(n, k);
-    for (size_t i=1; i <= fe.basis(n).size(); ++i)
-      EV[n](i) -= scale * conv * fe.basis(n)(i) * fe.detJxW;
+    EV[n].add(fe.basis(n), -scale*conv*fe.detJxW);
   }
 }