diff --git a/Apps/Common/SIMSolverAdap.h b/Apps/Common/SIMSolverAdap.h
index b38680b9..69152a6a 100644
--- a/Apps/Common/SIMSolverAdap.h
+++ b/Apps/Common/SIMSolverAdap.h
@@ -28,7 +28,11 @@ template<class T1> class SIMSolverAdap : public SIMSolver<T1>
 {
 public:
   //! \brief The constructor forwards to the parent class constructor.
-  SIMSolverAdap(T1& s1) : SIMSolver<T1>(s1), aSim(s1,false) {}
+  SIMSolverAdap(T1& s1) : SIMSolver<T1>(s1), aSim(s1,false)
+  {
+    this->S1.setSol(&aSim.getSolution());
+  }
+
   //! \brief Empty destructor.
   virtual ~SIMSolverAdap() {}
 
@@ -39,7 +43,6 @@ public:
     if (exporter)
       exporter->setNormPrefixes(aSim.getNormPrefixes());
 
-    this->S1.setSol(&aSim.getSolution());
     aSim.setupProjections();
     aSim.initAdaptor(0,2);
 
diff --git a/src/ASM/LR/ASMu2Drecovery.C b/src/ASM/LR/ASMu2Drecovery.C
index 1388c1f0..0eccf95d 100644
--- a/src/ASM/LR/ASMu2Drecovery.C
+++ b/src/ASM/LR/ASMu2Drecovery.C
@@ -414,19 +414,23 @@ LR::LRSplineSurface* ASMu2D::regularInterpolation (const RealArray& upar,
     return nullptr;
   }
 
-  DenseMatrix    A(nBasis,nBasis);
-  Matrix         B(points,true);
+  SparseMatrix   A(SparseMatrix::SUPERLU);
+  A.resize(nBasis, nBasis);
+  Matrix B2(points,true); // transpose to get one vector per field
+  StdVector B(B2);
   Go::BasisPtsSf splineValues;
 
   // Evaluate all basis functions at all points, stored in the A-matrix
   // (same row = same evaluation point)
   for (size_t i = 0; i < nBasis; i++)
   {
-    lrspline->computeBasis(upar[i],vpar[i],splineValues);
-    // optimization note:
-    // without an element id, splineValues will be stored as a full dense vector
-    for (size_t j = 0; j < nBasis; j++)
-      A(i+1,j+1) = splineValues.basisValues[j];
+    int iel = lrspline->getElementContaining(upar[i], vpar[i]);
+    lrspline->computeBasis(upar[i],vpar[i],splineValues, iel);
+    size_t k = 0;
+    for (auto& function : lrspline->getElement(iel)->support()) {
+      int j = function->getId();
+      A(i+1,j+1) = splineValues.basisValues[k++];
+    }
   }
 
   // Solve for all solution components - one right-hand-side for each
@@ -435,14 +439,17 @@ LR::LRSplineSurface* ASMu2D::regularInterpolation (const RealArray& upar,
 
   // Copy all basis functions and mesh
   LR::LRSplineSurface* ans = lrspline->copy();
-  ans->rebuildDimension(B.cols());
+  ans->rebuildDimension(points.rows());
 
-  // Swap around the control point values
-  for (size_t j = 0; j < B.cols(); j++) {
-    size_t i = 0;
-    for (LR::Basisfunction* b : ans->getAllBasisfunctions())
-      b->cp()[j] = B(++i,j+1);
+  // Back to interleaved data
+  std::vector<double> interleave;
+  interleave.reserve(B.dim());
+  for (size_t i = 0; i < nBasis; ++i)
+    for (size_t j = 0; j < points.rows(); j++) {
+        interleave.push_back(B(1+j*points.cols()+i));
   }
 
+  ans->setControlPoints(interleave);
+
   return ans;
 }