diff --git a/src/ASM/ASMmxBase.C b/src/ASM/ASMmxBase.C
index 5581d2eb..ce73847b 100644
--- a/src/ASM/ASMmxBase.C
+++ b/src/ASM/ASMmxBase.C
@@ -18,93 +18,88 @@
 #include "GoTools/trivariate/VolumeInterpolator.h"
 #include <numeric>
 
-char ASMmxBase::geoBasis            = 2;
+char ASMmxBase::geoBasis             = 2;
 ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
 
 
-ASMmxBase::ASMmxBase (const std::vector<unsigned char>& n_f)
-{
-  nfx = n_f;
-}
-
-
 void ASMmxBase::initMx (const std::vector<int>& MLGN, const int* sysMadof)
 {
-  MADOF.resize(MLGN.size());
-  for (size_t i = 0; i < MADOF.size(); i++)
-    MADOF[i] = sysMadof[MLGN[i]-1]-1;
+  MADOF.clear();
+  MADOF.reserve(MLGN.size());
+  for (int n : MLGN)
+    MADOF.push_back(sysMadof[n-1]-1);
 }
 
 
-void ASMmxBase::extractNodeVecMx (const Vector& globRes, Vector& nodeVec,
-				  int basis) const
+void ASMmxBase::extractNodeVecMx (const Vector& glbVec, Vector& nodVec,
+                                  int basis) const
 {
-  if (basis > (int)nfx.size())
-    basis = 0;
-
-  size_t len=0;
-  if (basis == 0)
-    for (size_t i=0;i<nfx.size();++i)
+  int b0 = basis-1, b1 = basis;
+  size_t ofs = 0, len = 0;
+  if (basis < 1 || basis > (int)nfx.size())
+  {
+    b0 = 0;
+    b1 = nfx.size();
+    for (int i = b0; i < b1; i++)
       len += nfx[i]*nb[i];
-  else
-    len = nfx[basis-1]*nb[basis-1];
-
-  nodeVec.resize(len);
-
-  size_t i, j;
-  int idof, ldof = 0;
-  size_t k=basis==0?1:basis;
-  size_t ofs = std::accumulate(nb.begin(), nb.begin()+k-1, 0);
-  for (; k < (basis==0?nfx.size()+1:(size_t)basis+1); ++k) {
-    for (i = ofs; i < nb[k-1]+ofs; i++)
-    {
-      idof = MADOF[i];
-      for (j = 0; j < nfx[k-1]; j++, ldof++)
-	nodeVec[ldof] = globRes[idof++];
-    }
-    ofs += nb[k-1];
   }
+  else
+  {
+    len = nfx[basis-1]*nb[basis-1];
+    if (basis > 1)
+      ofs = std::accumulate(nb.begin(), nb.begin()+basis-1, 0);
+  }
+
+  nodVec.resize(len);
+
+  int ldof = 0;
+  for (int b = b0; b < b1; b++, ofs += nb[b])
+    for (size_t i = ofs; i < nb[b]+ofs; i++)
+    {
+      int idof = MADOF[i];
+      for (size_t j = 0; j < nfx[b]; j++)
+        nodVec[ldof++] = glbVec[idof++];
+    }
 }
 
 
-void ASMmxBase::injectNodeVecMx (Vector& globRes, const Vector& nodeVec,
+void ASMmxBase::injectNodeVecMx (Vector& glbVec, const Vector& nodVec,
                                  int basis) const
 {
-  if (basis > (int)nfx.size())
-    basis = 0;
-
-  size_t i, j;
-  int idof, ldof = 0;
-  size_t k=basis==0?1:basis;
-  size_t ofs = std::accumulate(nb.begin(), nb.begin()+k-1, 0);
-  for (; k < (basis==0?nfx.size()+1:(size_t)basis+1); ++k) {
-    for (i = ofs; i < nb[k-1]+ofs; i++)
-    {
-      idof = MADOF[i];
-      for (j = 0; j < nfx[k-1]; j++, ldof++)
-	globRes[idof++] = nodeVec[ldof];
-    }
+  int b0 = basis-1, b1 = basis;
+  size_t ofs = 0;
+  if (basis < 1 || basis > (int)nfx.size())
+  {
+    b0 = 0;
+    b1 = nfx.size();
   }
+  else if (basis > 1)
+    ofs = std::accumulate(nb.begin(), nb.begin()+basis-1, 0);
+
+  int ldof = 0;
+  for (int b = b0; b < b1; b++)
+    for (size_t i = ofs; i < nb[b]+ofs; i++)
+    {
+      int idof = MADOF[i];
+      for (size_t j = 0; j < nfx[b]; j++)
+        glbVec[idof++] = nodVec[ldof++];
+    }
 }
 
 
 bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol,
-			       const std::vector<int>& nodes) const
+                               const std::vector<int>& nodes) const
 {
   if (nodes.empty()) return true;
 
-  int low, high, nvar;
-  if ((size_t)nodes.front() <= nb[0])
-  {
+  int nvar, low = 1, high = nb.front();
+  if (nodes.front() <= high)
     nvar = nfx[0];
-    low  = 1;
-    high = nb[0];
-  }
   else
   {
     nvar = nfx[1];
-    low  = nb[0]+1;
-    high = nb[0]+nb[1];
+    low += nb[0];
+    high += nb[1];
   }
 
   sField.resize(nvar,nodes.size());
@@ -112,15 +107,15 @@ bool ASMmxBase::getSolutionMx (Matrix& sField, const Vector& locSol,
     if (nodes[i] < low || nodes[i] > high)
     {
       std::cerr <<" *** ASMmxBase::getSolutionMx: Node #"<< nodes[i]
-		<<" is out of range ["<< low <<","<< high <<"]."<< std::endl;
+                <<" is out of range ["<< low <<","<< high <<"]."<< std::endl;
       return false;
     }
     else
     {
       int idof = nvar*(nodes[i]-1);
-      if (low > 1) idof += nfx[0]*nb[0];
+      if (low > 1) idof += nfx.front()*nb.front();
       for (int j = 0; j < nvar; j++)
-	sField(j+1,i+1) = locSol[idof++];
+        sField(j+1,i+1) = locSol[idof++];
     }
 
   return true;
diff --git a/src/ASM/ASMmxBase.h b/src/ASM/ASMmxBase.h
index 5d5237c7..115a709a 100644
--- a/src/ASM/ASMmxBase.h
+++ b/src/ASM/ASMmxBase.h
@@ -32,7 +32,7 @@ class ASMmxBase
 protected:
   //! \brief The constructor sets the number of field variables.
   //! \param[in] n_f Number of nodal variables in each field
-  explicit ASMmxBase(const std::vector<unsigned char>& n_f);
+  explicit ASMmxBase(const std::vector<unsigned char>& n_f) : nfx(n_f) {}
 
   //! \brief Initializes the patch level MADOF array.
   //! \param[in] MLGN Matrix of local-to-global node numbers
@@ -40,25 +40,23 @@ protected:
   void initMx(const std::vector<int>& MLGN, const int* sysMadof);
 
   //! \brief Extracts nodal results for this patch from the global vector.
-  //! \param[in] globVec Global solution vector in DOF-order
-  //! \param[out] nodeVec Nodal result vector for this patch
+  //! \param[in] glbVec Global solution vector in DOF-order
+  //! \param[out] nodVec Nodal result vector for this patch
   //! \param[in] basis Which basis to extract the nodal values for
-  void extractNodeVecMx(const Vector& globVec, Vector& nodeVec,
-			int basis = 0) const;
+  void extractNodeVecMx(const Vector& glbVec, Vector& nodVec, int basis) const;
 
   //! \brief Injects nodal results for this patch into a global vector.
-  //! \param[out] globVec Global solution vector in DOF-order
-  //! \param[in] nodeVec Nodal result vector for this patch
+  //! \param[out] glbVec Global solution vector in DOF-order
+  //! \param[in] nodVec Nodal result vector for this patch
   //! \param[in] basis Which basis to inject the nodal values for
-  void injectNodeVecMx(Vector& globVec, const Vector& nodeVec,
-                       int basis = 0) const;
+  void injectNodeVecMx(Vector& glbVec, const Vector& nodVec, int basis) const;
 
   //! \brief Extracts the primary solution field at the specified nodes.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector local to current patch
   //! \param[in] nodes 1-based local node numbers to extract solution for
   bool getSolutionMx(Matrix& sField, const Vector& locSol,
-		     const std::vector<int>& nodes) const;
+                     const std::vector<int>& nodes) const;
 
 public:
   //! \brief Enum defining available mixed formulation types.
@@ -95,7 +93,11 @@ private:
   std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch
 
 protected:
-  std::vector<size_t>        nb;  //!< Number of basis functions in each basis
+  //! Number of basis functions per element in each basis
+  std::vector<size_t> elem_size;
+  //! Total number of basis functions in each basis
+  std::vector<size_t> nb;
+
   std::vector<unsigned char> nfx; //!< Number of fields on each basis
 };
 
diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C
index df26084a..8af5946c 100644
--- a/src/ASM/ASMs2Dmx.C
+++ b/src/ASM/ASMs2Dmx.C
@@ -35,15 +35,13 @@
 #endif
 
 
-ASMs2Dmx::ASMs2Dmx (unsigned char n_s,
-		    const std::vector<unsigned char>& n_f)
+ASMs2Dmx::ASMs2Dmx (unsigned char n_s, const CharVec& n_f)
   : ASMs2D(n_s, *std::max_element(n_f.begin(),n_f.end())), ASMmxBase(n_f)
 {
 }
 
 
-ASMs2Dmx::ASMs2Dmx (const ASMs2Dmx& patch,
-                    const std::vector<unsigned char>& n_f)
+ASMs2Dmx::ASMs2Dmx (const ASMs2Dmx& patch, const CharVec& n_f)
   : ASMs2D(patch), ASMmxBase(n_f)
 {
   nb = patch.nb;
@@ -191,8 +189,13 @@ bool ASMs2Dmx::generateFEMTopology ()
   geo = surf = m_basis[geoBasis-1]->clone();
 
   nb.clear();
-  for (auto it : m_basis)
+  nb.reserve(m_basis.size());
+  elem_size.clear();
+  elem_size.reserve(m_basis.size());
+  for (auto& it : m_basis) {
     nb.push_back(it->numCoefs_u()*it->numCoefs_v());
+    elem_size.push_back(it->order_u()*it->order_v());
+  }
 
   if (!nodeInd.empty() && !shareFE)
   {
@@ -440,12 +443,8 @@ double ASMs2Dmx::getParametricArea (int iel) const
   if (MNPC[iel-1].empty())
     return 0.0;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
-
-  int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
-                                          elem_sizes.begin()+geoBasis, -1)];
+  int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
+                                          elem_size.begin()+geoBasis, -1)];
 #ifdef INDEX_CHECK
   if (inod1 < 0 || (size_t)inod1 >= nnod)
   {
@@ -475,12 +474,8 @@ double ASMs2Dmx::getParametricLength (int iel, int dir) const
   if (MNPC[iel-1].empty())
     return 0.0;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
-
-  int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
-                                          elem_sizes.begin()+geoBasis, -1)];
+  int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
+                                          elem_size.begin()+geoBasis, -1)];
 #ifdef INDEX_CHECK
   if (inod1 < 0 || (size_t)inod1 >= nnod)
   {
@@ -537,9 +532,6 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
   const int n1 = surf->numCoefs_u();
   const int nel1 = n1 - p1 + 1;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
 
   // === Assembly loop over all elements in the patch ==========================
 
@@ -547,7 +539,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
   for (size_t g=0;g<threadGroups.size() && ok;++g) {
 #pragma omp parallel for schedule(static)
     for (size_t t=0;t<threadGroups[g].size();++t) {
-      MxFiniteElement fe(elem_sizes);
+      MxFiniteElement fe(elem_size);
       std::vector<Matrix> dNxdu(m_basis.size());
       std::vector<Matrix3D> d2Nxdu2(m_basis.size());
       Matrix3D Hess;
@@ -590,8 +582,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
         }
 
         // Initialize element quantities
-        LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false);
-        if (!integrand.initElement(MNPC[iel-1],elem_sizes,nb,*A))
+        LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,false);
+        if (!integrand.initElement(MNPC[iel-1],elem_size,nb,*A))
         {
           A->destruct();
           ok = false;
@@ -726,19 +718,15 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
   const int n1 = surf->numCoefs_u();
   const int n2 = surf->numCoefs_v();
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
-
   std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
   size_t firstp = iit == firstBp.end() ? 0 : iit->second;
 
-  MxFiniteElement fe(elem_sizes);
+  MxFiniteElement fe(elem_size);
   fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
   fe.u = gpar[0](1,1);
   fe.v = gpar[1](1,1);
 
-  std::vector<Matrix> dNxdu(m_basis.size());
+  Matrices dNxdu(m_basis.size());
   Matrix Xnod, Jac;
   Vec4   X;
   Vec3   normal;
@@ -775,8 +763,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
         fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
 
       // Initialize element quantities
-      LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true);
-      bool ok = integrand.initElementBou(MNPC[iel-1],elem_sizes,nb,*A);
+      LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,true);
+      bool ok = integrand.initElementBou(MNPC[iel-1],elem_size,nb,*A);
 
 
       // --- Integration loop over all Gauss points along the edge -------------
@@ -859,11 +847,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
   const int n1 = surf->numCoefs_u();
   const int n2 = surf->numCoefs_v();
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
-  std::vector<size_t> elem_sizes2(elem_sizes);
-  std::copy(elem_sizes.begin(), elem_sizes.end(), std::back_inserter(elem_sizes2));
+  std::vector<size_t> elem_sizes2(elem_size);
+  std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
 
   MxFiniteElement fe(elem_sizes2);
   Matrix        dNdu, Xnod, Jac;
@@ -903,8 +888,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
         fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
 
       // Initialize element quantities
-      LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel);
-      bool ok = integrand.initElement(MNPC[iel],elem_sizes,nb,*A);
+      LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel);
+      bool ok = integrand.initElement(MNPC[iel],elem_size,nb,*A);
       size_t origSize = A->vec.size();
 
       // Loop over the element edges with contributions
@@ -924,8 +909,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
             kel += iedge > 3 ? n1-p1+1 : -(n1-p1+1);
 
           // initialize neighbor element
-          LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes,kel+1);
-          ok &= integrand.initElement(MNPC[kel],elem_sizes,nb,*A_neigh);
+          LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_size,kel+1);
+          ok &= integrand.initElement(MNPC[kel],elem_size,nb,*A_neigh);
           if (!A_neigh->vec.empty()) {
             A->vec.resize(origSize+A_neigh->vec.size());
             std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize);
@@ -961,7 +946,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
             }
 
             // Fetch basis function derivatives at current integration point
-            std::vector<Matrix> dNxdu(m_basis.size()*2);
+            Matrices dNxdu(m_basis.size()*2);
             for (size_t b = 0; b < m_basis.size(); ++b) {
               Go::BasisDerivsSf spline;
               this->getBasis(b+1)->computeBasis(fe.u, fe.v, spline, edgeDir < 0);
@@ -1120,17 +1105,13 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     }
   }
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order_u()*it->order_v());
-
   // Fetch nodal (control point) coordinates
   Matrix Xnod, Xtmp;
   this->getNodalCoordinates(Xnod);
 
-  MxFiniteElement fe(elem_sizes,firstIp);
+  MxFiniteElement fe(elem_size,firstIp);
   Vector          solPt;
-  std::vector<Matrix> dNxdu(m_basis.size());
+  Matrices        dNxdu(m_basis.size());
   Matrix          Jac;
   Vec3            X;
 
@@ -1139,7 +1120,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
   for (size_t i = 0; i < nPoints; i++, fe.iGP++)
   {
     // Fetch indices of the non-zero basis functions at this point
-    std::vector<IntVec> ip(m_basis.size());
+    IntMat ip(m_basis.size());
     IntVec ipa;
     size_t ofs = 0;
     for (size_t b = 0; b < m_basis.size(); ++b) {
@@ -1180,7 +1161,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     X = Xtmp * fe.basis(geoBasis);
 
     // Now evaluate the solution field
-    if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb))
+    if (!integrand.evalSol(solPt,fe,X,ipa,elem_size,nb))
       return false;
     else if (sField.empty())
       sField.resize(solPt.size(),nPoints,true);
diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C
index 2cd7cf3e..47a4b47d 100644
--- a/src/ASM/ASMs3Dmx.C
+++ b/src/ASM/ASMs3Dmx.C
@@ -35,14 +35,13 @@
 #endif
 
 
-ASMs3Dmx::ASMs3Dmx (const std::vector<unsigned char>& n_f)
+ASMs3Dmx::ASMs3Dmx (const CharVec& n_f)
   : ASMs3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f)
 {
 }
 
 
-ASMs3Dmx::ASMs3Dmx (const ASMs3Dmx& patch,
-                    const std::vector<unsigned char>& n_f)
+ASMs3Dmx::ASMs3Dmx (const ASMs3Dmx& patch, const CharVec& n_f)
   : ASMs3D(patch), ASMmxBase(n_f)
 {
   m_basis = patch.m_basis;
@@ -192,8 +191,13 @@ bool ASMs3Dmx::generateFEMTopology ()
   geo = svol = m_basis[geoBasis-1]->clone();
 
   nb.clear();
-  for (auto it : m_basis)
+  elem_size.clear();
+  nb.reserve(m_basis.size());
+  elem_size.reserve(m_basis.size());
+  for (auto& it : m_basis) {
     nb.push_back(it->numCoefs(0)*it->numCoefs(1)*it->numCoefs(2));
+    elem_size.push_back(it->order(0)*it->order(1)*it->order(2));
+  }
 
   if (!nodeInd.empty() && !shareFE)
   {
@@ -474,10 +478,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
     for (size_t i=0;i<m_basis.size();++i)
       m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]);
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-
   const int p1 = svol->order(0);
   const int p2 = svol->order(1);
   const int p3 = svol->order(2);
@@ -493,7 +493,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
   for (size_t g=0;g<threadGroupsVol.size() && ok;++g) {
 #pragma omp parallel for schedule(static)
     for (size_t t=0;t<threadGroupsVol[g].size();++t) {
-      MxFiniteElement fe(elem_sizes);
+      MxFiniteElement fe(elem_size);
       std::vector<Matrix> dNxdu(m_basis.size());
       std::vector<Matrix3D> d2Nxdu2(m_basis.size());
       Matrix3D Hess;
@@ -538,8 +538,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
         }
 
         // Initialize element quantities
-        LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false);
-        if (!integrand.initElement(MNPC[iel-1],elem_sizes,nb,*A))
+        LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,false);
+        if (!integrand.initElement(MNPC[iel-1],elem_size,nb,*A))
         {
           A->destruct();
           ok = false;
@@ -690,10 +690,6 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
   const int p2 = svol->order(1);
   const int p3 = svol->order(2);
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-
   const int nel1 = n1 - p1 + 1;
   const int nel2 = n2 - p2 + 1;
 
@@ -707,15 +703,15 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
   for (size_t g = 0; g < threadGrp.size() && ok; ++g) {
 #pragma omp parallel for schedule(static)
     for (size_t t = 0; t < threadGrp[g].size(); ++t) {
-      MxFiniteElement fe(elem_sizes);
+      MxFiniteElement fe(elem_size);
       fe.xi = fe.eta = fe.zeta = faceDir < 0 ? -1.0 : 1.0;
       fe.u = gpar[0](1,1);
       fe.v = gpar[1](1,1);
       fe.w = gpar[2](1,1);
 
-      std::vector<Matrix> dNxdu(m_basis.size());
+      Matrices dNxdu(m_basis.size());
       Matrix Xnod, Jac;
-      double   param[3] = { fe.u, fe.v, fe.w };
+      double param[3] = { fe.u, fe.v, fe.w };
       Vec4   X;
       Vec3   normal;
       for (size_t l = 0; l < threadGrp[g][t].size() && ok; ++l)
@@ -747,8 +743,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
           fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
 
 	// Initialize element quantities
-        LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true);
-	if (!integrand.initElementBou(MNPC[iel-1],elem_sizes,nb,*A))
+        LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel,true);
+	if (!integrand.initElementBou(MNPC[iel-1],elem_size,nb,*A))
         {
           A->destruct();
           ok = false;
@@ -867,11 +863,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
   const int nel1 = n1 - p1 + 1;
   const int nel2 = n2 - p2 + 1;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-  std::vector<size_t> elem_sizes2(elem_sizes);
-  std::copy(elem_sizes.begin(), elem_sizes.end(), std::back_inserter(elem_sizes2));
+  std::vector<size_t> elem_sizes2(elem_size);
+  std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
 
   MxFiniteElement fe(elem_sizes2);
   Matrix        dNdu, Xnod, Jac;
@@ -912,8 +905,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
           fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
 
         // Initialize element quantities
-        LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel);
-        bool ok = integrand.initElement(MNPC[iel],elem_sizes,nb,*A);
+        LocalIntegral* A = integrand.getLocalIntegral(elem_size,fe.iel);
+        bool ok = integrand.initElement(MNPC[iel],elem_size,nb,*A);
         size_t origSize = A->vec.size();
 
         // Loop over the element edges with contributions
@@ -943,8 +936,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
               kel += iface == 6 ? (n2-p2+1)*(n1-p1+1) : -(n2-p2+1)*(n1-p1+1);
 
             // initialize neighbor element
-            LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes,kel+1);
-            ok &= integrand.initElement(MNPC[kel],elem_sizes,nb,*A_neigh);
+            LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_size,kel+1);
+            ok &= integrand.initElement(MNPC[kel],elem_size,nb,*A_neigh);
             if (!A_neigh->vec.empty()) {
               A->vec.resize(origSize+A_neigh->vec.size());
               std::copy(A_neigh->vec.begin(), A_neigh->vec.end(), A->vec.begin()+origSize);
@@ -1007,7 +1000,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
                 }
 
                 // Fetch basis function derivatives at current integration point
-                std::vector<Matrix> dNxdu(m_basis.size()*2);
+                Matrices dNxdu(m_basis.size()*2);
                 for (size_t b = 0; b < m_basis.size(); ++b) {
                   Go::BasisDerivs spline;
                   this->getBasis(b+1)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0);
@@ -1172,17 +1165,13 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     }
   }
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-
   // Fetch nodal (control point) coordinates
   Matrix Xnod, Xtmp;
   this->getNodalCoordinates(Xnod);
 
-  MxFiniteElement fe(elem_sizes,firstIp);
+  MxFiniteElement fe(elem_size,firstIp);
   Vector          solPt;
-  std::vector<Matrix> dNxdu(m_basis.size());
+  Matrices        dNxdu(m_basis.size());
   Matrix          Jac;
   Vec3            X;
 
@@ -1191,7 +1180,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
   for (size_t i = 0; i < nPoints; i++, fe.iGP++)
   {
     // Fetch indices of the non-zero basis functions at this point
-    std::vector<IntVec> ip(m_basis.size());
+    IntMat ip(m_basis.size());
     IntVec ipa;
     size_t ofs = 0;
     for (size_t b = 0; b < m_basis.size(); ++b) {
@@ -1234,7 +1223,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     X = Xtmp * fe.basis(geoBasis);
 
     // Now evaluate the solution field
-    if (!integrand.evalSol(solPt,fe,X,ipa,elem_sizes,nb))
+    if (!integrand.evalSol(solPt,fe,X,ipa,elem_size,nb))
       return false;
     else if (sField.empty())
       sField.resize(solPt.size(),nPoints,true);
@@ -1283,12 +1272,8 @@ double ASMs3Dmx::getParametricVolume (int iel) const
   if (MNPC[iel-1].empty())
     return 0.0;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-
-  int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
-                                          elem_sizes.begin()+geoBasis, -1)];
+  int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
+                                          elem_size.begin()+geoBasis, -1)];
 #ifdef INDEX_CHECK
   if (inod1 < 0 || (size_t)inod1 >= nnod)
   {
@@ -1318,12 +1303,8 @@ double ASMs3Dmx::getParametricArea (int iel, int dir) const
   if (MNPC[iel-1].empty())
     return 0.0;
 
-  std::vector<size_t> elem_sizes;
-  for (auto& it : m_basis)
-    elem_sizes.push_back(it->order(0)*it->order(1)*it->order(2));
-
-  int inod1 = MNPC[iel-1][std::accumulate(elem_sizes.begin(),
-                                          elem_sizes.begin()+geoBasis, -1)];
+  int inod1 = MNPC[iel-1][std::accumulate(elem_size.begin(),
+                                          elem_size.begin()+geoBasis, -1)];
 #ifdef INDEX_CHECK
   if (inod1 < 0 || (size_t)inod1 >= nnod)
   {