diff --git a/configure b/configure
index e28176191bd..78f77119e6f 100755
--- a/configure
+++ b/configure
@@ -63450,13 +63450,13 @@ ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
           cat confdefs.h - <<_ACEOF >conftest.$ac_ext
 /* end confdefs.h.  */
 
-          #include "metaphysicl/numberarray.h"
+          #include "metaphysicl/raw_type.h"
 
 int
 main (void)
 {
 
-              MetaPhysicL::NumberArray<4, double> x;
+              MetaPhysicL::RawType<double> x;
 
   ;
   return 0;
diff --git a/contrib/fparser/fparser_ad.cc b/contrib/fparser/fparser_ad.cc
index 1bfe2c10738..2720e0cc5ac 100644
--- a/contrib/fparser/fparser_ad.cc
+++ b/contrib/fparser/fparser_ad.cc
@@ -800,77 +800,77 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
       case cDiv:
         --sp; ccout << "s[" << sp << "] /= s[" << (sp+1) << "];\n"; break;
       case cMod:
-        --sp; ccout << "s[" << sp << "] = std::fmod(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
+        --sp; ccout << "using std::fmod; s[" << sp << "] = fmod(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
       case cRDiv:
         --sp; ccout << "s[" << sp << "] = s[" << (sp+1) << "] / s[" << sp << "];\n"; break;
 
       case cSin:
-        ccout << "s[" << sp << "] = std::sin(s[" << sp << "]);\n"; break;
+        ccout << "using std::sin; s[" << sp << "] = sin(s[" << sp << "]);\n"; break;
       case cCos:
-        ccout << "s[" << sp << "] = std::cos(s[" << sp << "]);\n"; break;
+        ccout << "using std::cos; s[" << sp << "] = cos(s[" << sp << "]);\n"; break;
       case cTan:
-        ccout << "s[" << sp << "] = std::tan(s[" << sp << "]);\n"; break;
+        ccout << "using std::tan; s[" << sp << "] = tan(s[" << sp << "]);\n"; break;
       case cSinh:
-        ccout << "s[" << sp << "] = std::sinh(s[" << sp << "]);\n"; break;
+        ccout << "using std::sinh; s[" << sp << "] = sinh(s[" << sp << "]);\n"; break;
       case cCosh:
-        ccout << "s[" << sp << "] = std::cosh(s[" << sp << "]);\n"; break;
+        ccout << "using std::cosh; s[" << sp << "] = cosh(s[" << sp << "]);\n"; break;
       case cTanh:
-        ccout << "s[" << sp << "] = std::tanh(s[" << sp << "]);\n"; break;
+        ccout << "using std::tanh; s[" << sp << "] = tanh(s[" << sp << "]);\n"; break;
       // TODO: div by zero -> this->mData->mEvalErrorType=1; return Value_t(0);
       case cCsc:
-        ccout << "s[" << sp << "] = 1.0/std::sin(s[" << sp << "]);\n"; break;
+        ccout << "using std::sin; s[" << sp << "] = 1.0/sin(s[" << sp << "]);\n"; break;
       case cSec:
-        ccout << "s[" << sp << "] = 1.0/std::cos(s[" << sp << "]);\n"; break;
+        ccout << "using std::cos; s[" << sp << "] = 1.0/cos(s[" << sp << "]);\n"; break;
       case cCot:
-        ccout << "s[" << sp << "] = 1.0/std::tan(s[" << sp << "]);\n"; break;
+        ccout << "using std::tan; s[" << sp << "] = 1.0/tan(s[" << sp << "]);\n"; break;
       case cSinCos:
-        ccout << "s[" << (sp+1) << "] = std::cos(s[" << sp << "]);\n";
-        ccout << "s[" << sp << "] = std::sin(s[" << sp << "]);\n";
+        ccout << "using std::cos; s[" << (sp+1) << "] = cos(s[" << sp << "]);\n";
+        ccout << "using std::sin; s[" << sp << "] = sin(s[" << sp << "]);\n";
         ++sp;
         break;
       case cSinhCosh:
-        ccout << "s[" << (sp+1) << "] = std::cosh(s[" << sp << "]);\n";
-        ccout << "s[" << sp << "] = std::sinh(s[" << sp << "]);\n";
+        ccout << "using std::cosh; s[" << (sp+1) << "] = cosh(s[" << sp << "]);\n";
+        ccout << "using std::sinh; s[" << sp << "] = sinh(s[" << sp << "]);\n";
         ++sp;
         break;
       case cAsin:
-        ccout << "s[" << sp << "] = std::asin(s[" << sp << "]);\n"; break;
+        ccout << " using std::asin; s[" << sp << "] = asin(s[" << sp << "]);\n"; break;
       case cAcos:
-        ccout << "s[" << sp << "] = std::acos(s[" << sp << "]);\n"; break;
+        ccout << " using std::acos; s[" << sp << "] = acos(s[" << sp << "]);\n"; break;
       case cAsinh:
-        ccout << "s[" << sp << "] = std::asinh(s[" << sp << "]);\n"; break;
+        ccout << " using std::asinh; s[" << sp << "] = asinh(s[" << sp << "]);\n"; break;
       case cAcosh:
-        ccout << "s[" << sp << "] = std::acosh(s[" << sp << "]);\n"; break;
+        ccout << " using std::acosh; s[" << sp << "] = acosh(s[" << sp << "]);\n"; break;
       case cAtan:
-        ccout << "s[" << sp << "] = std::atan(s[" << sp << "]);\n"; break;
+        ccout << " using std::atan; s[" << sp << "] = atan(s[" << sp << "]);\n"; break;
       case cAtanh:
-        ccout << "s[" << sp << "] = std::atanh(s[" << sp << "]);\n"; break;
+        ccout << " using std::atanh; s[" << sp << "] = atanh(s[" << sp << "]);\n"; break;
       case cAtan2:
-        --sp; ccout << "s[" << sp << "] = std::atan2(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
+        --sp; ccout << " using std::atan2; s[" << sp << "] = atan2(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
       case cHypot:
-        --sp; ccout << "s[" << sp << "] = std::sqrt(s[" << sp << "]*s[" << sp << "] + s[" << (sp+1) << "]*s[" << (sp+1) << "]);\n"; break;
+        --sp; ccout << " using std::sqrt; s[" << sp << "] = sqrt(s[" << sp << "]*s[" << sp << "] + s[" << (sp+1) << "]*s[" << (sp+1) << "]);\n"; break;
 
       case cAbs:
-        ccout << "s[" << sp << "] = std::abs(s[" << sp << "]);\n"; break;
+        ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]);\n"; break;
       case cMax:
         --sp; ccout << "s[" << sp << "] = s[" << sp << "] > s[" << (sp+1) << "] ? s[" << sp << "] : s[" << (sp+1) << "];\n"; break;
       case cMin:
         --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? s[" << sp << "] : s[" << (sp+1) << "];\n"; break;
       case cTrunc:
-        ccout << "s[" << sp << "] = s[" << sp << "] < 0 ? std::ceil(s[" << sp << "]) : std::floor(s[" << sp << "]);\n"; break;
+        ccout << "using std::ceil, std::floor; s[" << sp << "] = s[" << sp << "] < 0 ? ceil(s[" << sp << "]) : floor(s[" << sp << "]);\n"; break;
       case cCeil:
-        ccout << "s[" << sp << "] = std::ceil(s[" << sp << "]);\n"; break;
+        ccout << "using std::ceil; s[" << sp << "] = ceil(s[" << sp << "]);\n"; break;
       case cFloor:
-        ccout << "s[" << sp << "] = std::floor(s[" << sp << "]);\n"; break;
+        ccout << "using std::floor; s[" << sp << "] = floor(s[" << sp << "]);\n"; break;
       case cInt:
-        ccout << "s[" << sp << "] = s[" << sp << "] < 0 ? std::ceil(s[" << sp << "] - 0.5) : std::floor(s[" << sp << "] + 0.5);\n"; break;
+        ccout << "using std::ceil, std::floor; s[" << sp << "] = s[" << sp << "] < 0 ? ceil(s[" << sp << "] - 0.5) : floor(s[" << sp << "] + 0.5);\n"; break;
 
       case cEqual:
         //--sp; ccout << "s[" << sp << "] = s[" << sp << "] == s[" << (sp+1) << "];\n"; break;
-        --sp; ccout << "s[" << sp << "] = std::abs(s[" << sp << "] - s[" << (sp+1) << "]) <= eps;\n"; break;
+        --sp; ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "] - s[" << (sp+1) << "]) <= eps;\n"; break;
       case cNEqual:
         //--sp; ccout << "s[" << sp << "] = s[" << sp << "] != s[" << (sp+1) << "];\n"; break;
-        --sp; ccout << "s[" << sp << "] = std::abs(s[" << sp << "] - s[" << (sp+1) << "]) > eps;\n"; break;
+        --sp; ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "] - s[" << (sp+1) << "]) > eps;\n"; break;
       case cLess:
         --sp; ccout << "s[" << sp << "] = s[" << sp << "] < (s[" << (sp+1) << "] - eps);\n"; break;
       case cLessOrEq:
@@ -880,33 +880,33 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
       case cGreaterOrEq:
         --sp; ccout << "s[" << sp << "] = (s[" << sp << "] + eps) >= s[" << (sp+1) << "];\n"; break;
       case cNot:
-        ccout << "s[" << sp << "] = std::abs(s[" << sp << "]) < 0.5;\n"; break;
+        ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]) < 0.5;\n"; break;
       case cNotNot:
-        ccout << "s[" << sp << "] = std::abs(s[" << sp << "]) >= 0.5;\n"; break;
+        ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]) >= 0.5;\n"; break;
       case cAbsNot:
         ccout << "s[" << sp << "] = s[" << sp << "] < 0.5;\n"; break;
       case cAbsNotNot:
         ccout << "s[" << sp << "] = s[" << sp << "] >= 0.5;\n"; break;
       case cOr:
-        --sp; ccout << "s[" << sp << "] = (std::abs(s[" << sp << "]) >= 0.5) || (std::abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
+        --sp; ccout << "using std::abs; s[" << sp << "] = (abs(s[" << sp << "]) >= 0.5) || (abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
       case cAbsOr:
         --sp; ccout << "s[" << sp << "] = (s[" << sp << "] >= 0.5) || (s[" << (sp+1) << "] >= 0.5);\n"; break;
       case cAnd:
-        --sp; ccout << "s[" << sp << "] = (std::abs(s[" << sp << "]) >= 0.5) && (std::abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
+        --sp; ccout << "using std::abs; s[" << sp << "] = (abs(s[" << sp << "]) >= 0.5) && (abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
       case cAbsAnd:
         --sp; ccout << "s[" << sp << "] = (s[" << sp << "] >= 0.5) && (s[" << (sp+1) << "] >= 0.5);\n"; break;
 
       case cLog:
-        ccout << "s[" << sp << "] = std::log(s[" << sp << "]);\n"; break;
+        ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "]);\n"; break;
       case cLog2:
 #ifdef FP_SUPPORT_CPLUSPLUS11_MATH_FUNCS
-        ccout << "s[" << (sp-1) << "] = std::log2(s[" << (sp-1) << "]);\n";
+        ccout << "using std::log2; s[" << (sp-1) << "] = log2(s[" << (sp-1) << "]);\n";
 #else
-        ccout << "s[" << sp << "] = std::log(s[" << sp << "])/log(2.0);\n";
+        ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "])/log(2.0);\n";
 #endif
         break;
       case cLog10:
-        ccout << "s[" << sp << "] = std::log10(s[" << sp << "]);\n"; break;
+        ccout << "using std::log10; s[" << sp << "] = log10(s[" << sp << "]);\n"; break;
 
       case cNeg:
         ccout << "s[" << sp << "] = -s[" << sp << "];\n"; break;
@@ -929,14 +929,14 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
         {
           // --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "]) + (s[" << sp << "] - s[" << (sp+1) << "]) / s[" << (sp+1) << "] : std::log(s[" << sp << "]);\n";
           // --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "]) - 1.5 + 2.0/s[" << (sp+1) << "] * s[" << sp << "] - 0.5/(s[" << (sp+1) << "]*s[" << (sp+1) << "]) * s[" << sp << "]*s[" << sp << "] : std::log(s[" << sp << "]);\n";
-          --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "])  +  (s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "]  -  std::pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],2.0)/2.0  +  std::pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],3.0)/3.0 : std::log(s[" << sp << "]);\n";
+          --sp; ccout << "using std::log, std::pow; s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? log(s[" << (sp+1) << "])  +  (s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "]  -  pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],2.0)/2.0  +  pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],3.0)/3.0 : log(s[" << sp << "]);\n";
         }
         else if (function == mFErf)
         {
 #if LIBMESH_HAVE_CXX11_ERF
-          ccout << "s[" << sp << "] = std::erf(s[" << sp << "]);\n";
+          ccout << "using std::erf; s[" << sp << "] = erf(s[" << sp << "]);\n";
 #else
-          std::cerr << "Libmesh is not compiled with c++11 so std::erf is not supported by JIT.\n";
+          std::cerr << "Libmesh is not compiled with c++11 so erf is not supported by JIT.\n";
           return false;
 #endif
         }
@@ -958,7 +958,7 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
         break;
       }
       case cLog2by:
-        --sp; ccout << "s[" << sp << "] = std::log(s[" << sp << "])/log(2.0) * s[" << (sp+1) << "];\n"; break;
+        --sp; ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "])/log(2.0) * s[" << (sp+1) << "];\n"; break;
       case cNop:
         break;
 #endif
@@ -966,20 +966,20 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
       case cSqr:
         ccout << "s[" << sp << "] *= s[" << sp << "];\n"; break;
       case cSqrt:
-        ccout << "s[" << sp << "] = std::sqrt(s[" << sp << "]);\n"; break;
+        ccout << "using std::sqrt; s[" << sp << "] = sqrt(s[" << sp << "]);\n"; break;
       case cRSqrt:
-        ccout << "s[" << sp << "] = std::pow(s[" << sp << "], (-0.5));\n"; break;
+        ccout << "using std::pow; s[" << sp << "] = pow(s[" << sp << "], (-0.5));\n"; break;
       case cPow:
-        --sp; ccout << "s[" << sp << "] = std::pow(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
+        --sp; ccout << "using std::pow; s[" << sp << "] = pow(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
       case cExp:
-        ccout << "s[" << sp << "] = std::exp(s[" << sp << "]);\n"; break;
+        ccout << "using std::exp; s[" << sp << "] = exp(s[" << sp << "]);\n"; break;
       case cExp2:
-        ccout << "s[" << sp << "] = std::pow(2.0, s[" << sp << "]);\n"; break;
+        ccout << "using std::pow; s[" << sp << "] = pow(2.0, s[" << sp << "]);\n"; break;
       case cCbrt:
 #ifdef FP_SUPPORT_CPLUSPLUS11_MATH_FUNCS
-        ccout << "s[" << sp << "] = std::cbrt(s[" << sp << "]);\n"; break;
+        ccout << "using std::cbrt; s[" << sp << "] = cbrt(s[" << sp << "]);\n"; break;
 #else
-        ccout << "s[" << sp << "] = s[" << sp << "] == 0 ? 0 : (s[" << sp << "] > 0 ? std::exp(std::log(s[" << sp << "])/3.0) : -std::exp(std::log(-s[" << sp << "])/3.0));\n"; break;
+        ccout << "using std::exp, std::log; s[" << sp << "] = s[" << sp << "] == 0 ? 0 : (s[" << sp << "] > 0 ? exp(log(s[" << sp << "])/3.0) : -exp(log(-s[" << sp << "])/3.0));\n"; break;
 #endif
 
       case cJump:
@@ -991,7 +991,7 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
         if (op == cIf)
           ccout << "if (s[" << sp-- << "] < 0.5) ";
         if (op == cAbsIf)
-          ccout << "if (std::abs(s[" << sp-- << "]) < 0.5) ";
+          ccout << "using std::abs; if (abs(s[" << sp-- << "]) < 0.5) ";
 
         if (ip >= ByteCode.size())
           ccout << "*ret = s[" << sp << "]; return;\n";
diff --git a/contrib/metaphysicl b/contrib/metaphysicl
index 34f1376858e..32f9b320684 160000
--- a/contrib/metaphysicl
+++ b/contrib/metaphysicl
@@ -1 +1 @@
-Subproject commit 34f1376858e9bd51d7ba53b4c3e59502890d31f3
+Subproject commit 32f9b32068496e722c496744cf641a85680c9fcd
diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h
index bd27382e06d..c38e8b9bbac 100644
--- a/include/numerics/dense_matrix.h
+++ b/include/numerics/dense_matrix.h
@@ -324,7 +324,7 @@ class DenseMatrix : public DenseMatrixBase<T>
    * This is the natural matrix norm that is compatible to the l1-norm
    * for vectors, i.e. \f$ |Mv|_1 \leq |M|_1 |v|_1 \f$.
    */
-  auto l1_norm () const -> decltype(std::abs(T(0)));
+  auto l1_norm () const;
 
   /**
    * \returns The linfty-norm of the matrix, that is, the max row sum:
@@ -334,7 +334,7 @@ class DenseMatrix : public DenseMatrixBase<T>
    * This is the natural matrix norm that is compatible to the
    * linfty-norm of vectors, i.e. \f$ |Mv|_\infty \leq |M|_\infty |v|_\infty \f$.
    */
-  auto linfty_norm () const -> decltype(std::abs(T(0)));
+  auto linfty_norm () const;
 
   /**
    * Left multiplies by the transpose of the matrix \p A.
@@ -1122,15 +1122,16 @@ auto DenseMatrix<T>::max () const -> decltype(libmesh_real(T(0)))
 
 template<typename T>
 inline
-auto DenseMatrix<T>::l1_norm () const -> decltype(std::abs(T(0)))
+auto DenseMatrix<T>::l1_norm () const
 {
   libmesh_assert (this->_m);
   libmesh_assert (this->_n);
 
-  auto columnsum = std::abs(T(0));
+  using std::abs;
+  auto columnsum = abs(T(0));
   for (unsigned int i=0; i!=this->_m; i++)
     {
-      columnsum += std::abs((*this)(i,0));
+      columnsum += abs((*this)(i,0));
     }
   auto my_max = columnsum;
   for (unsigned int j=1; j!=this->_n; j++)
@@ -1138,7 +1139,7 @@ auto DenseMatrix<T>::l1_norm () const -> decltype(std::abs(T(0)))
       columnsum = 0.;
       for (unsigned int i=0; i!=this->_m; i++)
         {
-          columnsum += std::abs((*this)(i,j));
+          columnsum += abs((*this)(i,j));
         }
       my_max = (my_max > columnsum? my_max : columnsum);
     }
@@ -1149,15 +1150,16 @@ auto DenseMatrix<T>::l1_norm () const -> decltype(std::abs(T(0)))
 
 template<typename T>
 inline
-auto DenseMatrix<T>::linfty_norm () const -> decltype(std::abs(T(0)))
+auto DenseMatrix<T>::linfty_norm () const
 {
   libmesh_assert (this->_m);
   libmesh_assert (this->_n);
+  using std::abs;
 
-  auto rowsum = std::abs(T(0));
+  auto rowsum = abs(T(0));
   for (unsigned int j=0; j!=this->_n; j++)
     {
-      rowsum += std::abs((*this)(0,j));
+      rowsum += abs((*this)(0,j));
     }
   auto my_max = rowsum;
   for (unsigned int i=1; i!=this->_m; i++)
@@ -1165,7 +1167,7 @@ auto DenseMatrix<T>::linfty_norm () const -> decltype(std::abs(T(0)))
       rowsum = 0.;
       for (unsigned int j=0; j!=this->_n; j++)
         {
-          rowsum += std::abs((*this)(i,j));
+          rowsum += abs((*this)(i,j));
         }
       my_max = (my_max > rowsum? my_max : rowsum);
     }
diff --git a/include/numerics/dense_matrix_impl.h b/include/numerics/dense_matrix_impl.h
index adeedb0ec3d..ee6decbfd27 100644
--- a/include/numerics/dense_matrix_impl.h
+++ b/include/numerics/dense_matrix_impl.h
@@ -630,6 +630,8 @@ void DenseMatrix<T>::_lu_back_substitute (const DenseVector<T> & b,
 template<typename T>
 void DenseMatrix<T>::_lu_decompose ()
 {
+  using std::abs;
+
   // If this function was called, there better not be any
   // previous decomposition of the matrix.
   libmesh_assert_equal_to (this->_decomposition_type, NONE);
@@ -648,11 +650,11 @@ void DenseMatrix<T>::_lu_decompose ()
       // Find the pivot row by searching down the i'th column
       _pivots[i] = i;
 
-      // std::abs(complex) must return a Real!
-      auto the_max = std::abs( A(i,i) );
+      // abs(complex) must return a Real!
+      auto the_max = abs( A(i,i) );
       for (unsigned int j=i+1; j<n_rows; ++j)
         {
-          auto candidate_max = std::abs( A(j,i) );
+          auto candidate_max = abs( A(j,i) );
           if (the_max < candidate_max)
             {
               the_max = candidate_max;
@@ -872,6 +874,8 @@ void DenseMatrix<T>::cholesky_solve (const DenseVector<T2> & b,
 template<typename T>
 void DenseMatrix<T>::_cholesky_decompose ()
 {
+  using std::sqrt;
+
   // If we called this function, there better not be any
   // previous decomposition of the matrix.
   libmesh_assert_equal_to (this->_decomposition_type, NONE);
@@ -901,7 +905,7 @@ void DenseMatrix<T>::_cholesky_decompose ()
                                    "Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
 #endif
 
-              A(i,i) = std::sqrt(A(i,j));
+              A(i,i) = sqrt(A(i,j));
             }
           else
             A(j,i) = A(i,j) / A(i,i);
diff --git a/include/numerics/tensor_tools.h b/include/numerics/tensor_tools.h
index 239d4ec25e6..a92ad862c8b 100644
--- a/include/numerics/tensor_tools.h
+++ b/include/numerics/tensor_tools.h
@@ -71,61 +71,61 @@ inner_product(const TypeNTensor<N,T> & a, const TypeNTensor<N,T2> & b)
 
 template<typename T>
 inline
-auto norm(const T & a) -> decltype(std::abs(a))
-{ return std::abs(a); }
+auto norm(const T & a)
+{ using std::abs; return abs(a); }
 
 template<typename T>
 inline
-T norm(std::complex<T> a) { return std::abs(a); }
+T norm(std::complex<T> a) { using std::abs; return abs(a); }
 
 template <typename T>
 inline
 auto norm(const TypeVector<T> & a) -> decltype(TensorTools::norm(T()))
-{return std::sqrt(a.norm_sq());}
+{using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
 inline
 auto norm(const VectorValue<T> & a) -> decltype(TensorTools::norm(T()))
-{return std::sqrt(a.norm_sq());}
+{using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
 inline
 auto norm(const TypeTensor<T> & a) -> decltype(TensorTools::norm(T()))
-{return std::sqrt(a.norm_sq());}
+{using std::sqrt; return sqrt(a.norm_sq());}
 
 template <typename T>
 inline
 auto norm(const TensorValue<T> & a) -> decltype(TensorTools::norm(T()))
-{return std::sqrt(a.norm_sq());}
+{using std::sqrt; return sqrt(a.norm_sq());}
 
 
 template<typename T>
 inline
-auto norm_sq(const T & a) -> decltype(std::norm(a))
-{ return std::norm(a); }
+auto norm_sq(const T & a)
+{ using std::norm; return norm(a); }
 
 template<typename T>
 inline
-T norm_sq(std::complex<T> a) { return std::norm(a); }
+T norm_sq(std::complex<T> a) { using std::norm; return norm(a); }
 
 template <typename T>
 inline
-auto norm_sq(const TypeVector<T> & a) -> decltype(std::norm(T()))
+auto norm_sq(const TypeVector<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
 inline
-auto norm_sq(const VectorValue<T> & a) -> decltype(std::norm(T()))
+auto norm_sq(const VectorValue<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
 inline
-auto norm_sq(const TypeTensor<T> & a) -> decltype(std::norm(T()))
+auto norm_sq(const TypeTensor<T> & a)
 {return a.norm_sq();}
 
 template <typename T>
 inline
-auto norm_sq(const TensorValue<T> & a) -> decltype(std::norm(T()))
+auto norm_sq(const TensorValue<T> & a)
 {return a.norm_sq();}
 
 
diff --git a/include/numerics/type_n_tensor.h b/include/numerics/type_n_tensor.h
index 5f00088db98..15bebec4f49 100644
--- a/include/numerics/type_n_tensor.h
+++ b/include/numerics/type_n_tensor.h
@@ -218,7 +218,7 @@ class TypeNTensor
    * \returns The Frobenius norm of the tensor, i.e. the square-root of
    * the sum of the elements squared.
    */
-  auto norm() const -> decltype(std::norm(T()))
+  auto norm() const
   {
     libmesh_not_implemented();
     return 0.;
@@ -228,7 +228,7 @@ class TypeNTensor
    * \returns The Frobenius norm of the tensor squared, i.e. the sum of the
    * entry magnitudes squared.
    */
-  auto norm_sq() const -> decltype(std::norm(T()))
+  auto norm_sq() const
   {
     libmesh_not_implemented();
     return 0.;
diff --git a/include/numerics/type_tensor.h b/include/numerics/type_tensor.h
index be3c0661917..67a79b430c4 100644
--- a/include/numerics/type_tensor.h
+++ b/include/numerics/type_tensor.h
@@ -364,13 +364,13 @@ class TypeTensor
    * \returns The Frobenius norm of the tensor, i.e. the square-root of
    * the sum of the elements squared.
    */
-  auto norm() const -> decltype(std::norm(T()));
+  auto norm() const;
 
   /**
    * \returns The Frobenius norm of the tensor squared, i.e. sum of the
    * element magnitudes squared.
    */
-  auto norm_sq() const -> decltype(std::norm(T()));
+  auto norm_sq() const;
 
   /**
    * \returns True if all values in the tensor are zero
@@ -1277,9 +1277,10 @@ TypeTensor<T>::contract (const TypeTensor<T2> & t) const
 
 template <typename T>
 inline
-auto TypeTensor<T>::norm() const -> decltype(std::norm(T()))
+auto TypeTensor<T>::norm() const
 {
-  return std::sqrt(this->norm_sq());
+  using std::sqrt;
+  return sqrt(this->norm_sq());
 }
 
 
@@ -1345,7 +1346,7 @@ void TypeTensor<T>::zero()
 
 template <typename T>
 inline
-auto TypeTensor<T>::norm_sq () const -> decltype(std::norm(T()))
+auto TypeTensor<T>::norm_sq () const
 {
   Real sum = 0.;
   for (unsigned int i=0; i<LIBMESH_DIM*LIBMESH_DIM; i++)
diff --git a/include/numerics/type_vector.h b/include/numerics/type_vector.h
index f12085cd5b0..d7d8f4ce24e 100644
--- a/include/numerics/type_vector.h
+++ b/include/numerics/type_vector.h
@@ -298,13 +298,13 @@ class TypeVector
    * \returns The magnitude of the vector, i.e. the square-root of the
    * sum of the elements squared.
    */
-  auto norm() const -> decltype(std::norm(T()));
+  auto norm() const;
 
   /**
    * \returns The magnitude of the vector squared, i.e. the sum of the
    * element magnitudes squared.
    */
-  auto norm_sq() const -> decltype(std::norm(T()));
+  auto norm_sq() const;
 
   /**
    * \returns The L1 norm of the vector
@@ -904,9 +904,10 @@ TypeVector<T>::cross(const TypeVector<T2> & p) const
 
 template <typename T>
 inline
-auto TypeVector<T>::norm() const -> decltype(std::norm(T()))
+auto TypeVector<T>::norm() const
 {
-  return std::sqrt(this->norm_sq());
+  using std::sqrt;
+  return sqrt(this->norm_sq());
 }
 
 
@@ -923,7 +924,7 @@ void TypeVector<T>::zero()
 
 template <typename T>
 inline
-auto TypeVector<T>::norm_sq() const -> decltype(std::norm(T()))
+auto TypeVector<T>::norm_sq() const
 {
 #if LIBMESH_DIM == 1
   return (TensorTools::norm_sq(_coords[0]));
diff --git a/m4/libmesh_metaphysicl.m4 b/m4/libmesh_metaphysicl.m4
index d4b9bca4d69..934bceb605d 100644
--- a/m4/libmesh_metaphysicl.m4
+++ b/m4/libmesh_metaphysicl.m4
@@ -73,9 +73,9 @@ AC_DEFUN([CONFIGURE_METAPHYSICL],
           CXXFLAGS="$CXXFLAGS $METAPHYSICL_INC"
 
           AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
-          @%:@include "metaphysicl/numberarray.h"
+          @%:@include "metaphysicl/raw_type.h"
           ]], [[
-              MetaPhysicL::NumberArray<4, double> x;
+              MetaPhysicL::RawType<double> x;
           ]])],[
               enablemetaphysicl=yes
               AC_MSG_RESULT(yes)