From cda28086408fb87679f4bbdb1f39abbeee9cbaba Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Wed, 9 Nov 2016 14:49:47 +0100
Subject: [PATCH] fixed: do not use a dense matrix for interpolation in ASMu2D

this is a sparse system
---
 src/ASM/LR/ASMu2Drecovery.C | 33 ++++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 13 deletions(-)

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;
 }