Merge pull request #69 from akva2/fix_lr_interpolation_hdf5

Fix two LR issues
This commit is contained in:
Kjetil Andre Johannessen 2016-11-10 15:14:11 +01:00 committed by GitHub
commit 8d81496e68
2 changed files with 25 additions and 15 deletions

View File

@ -28,7 +28,11 @@ template<class T1> class SIMSolverAdap : public SIMSolver<T1>
{ {
public: public:
//! \brief The constructor forwards to the parent class constructor. //! \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. //! \brief Empty destructor.
virtual ~SIMSolverAdap() {} virtual ~SIMSolverAdap() {}
@ -39,7 +43,6 @@ public:
if (exporter) if (exporter)
exporter->setNormPrefixes(aSim.getNormPrefixes()); exporter->setNormPrefixes(aSim.getNormPrefixes());
this->S1.setSol(&aSim.getSolution());
aSim.setupProjections(); aSim.setupProjections();
aSim.initAdaptor(0,2); aSim.initAdaptor(0,2);

View File

@ -414,19 +414,23 @@ LR::LRSplineSurface* ASMu2D::regularInterpolation (const RealArray& upar,
return nullptr; return nullptr;
} }
DenseMatrix A(nBasis,nBasis); SparseMatrix A(SparseMatrix::SUPERLU);
Matrix B(points,true); A.resize(nBasis, nBasis);
Matrix B2(points,true); // transpose to get one vector per field
StdVector B(B2);
Go::BasisPtsSf splineValues; Go::BasisPtsSf splineValues;
// Evaluate all basis functions at all points, stored in the A-matrix // Evaluate all basis functions at all points, stored in the A-matrix
// (same row = same evaluation point) // (same row = same evaluation point)
for (size_t i = 0; i < nBasis; i++) for (size_t i = 0; i < nBasis; i++)
{ {
lrspline->computeBasis(upar[i],vpar[i],splineValues); int iel = lrspline->getElementContaining(upar[i], vpar[i]);
// optimization note: lrspline->computeBasis(upar[i],vpar[i],splineValues, iel);
// without an element id, splineValues will be stored as a full dense vector size_t k = 0;
for (size_t j = 0; j < nBasis; j++) for (auto& function : lrspline->getElement(iel)->support()) {
A(i+1,j+1) = splineValues.basisValues[j]; int j = function->getId();
A(i+1,j+1) = splineValues.basisValues[k++];
}
} }
// Solve for all solution components - one right-hand-side for each // 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 // Copy all basis functions and mesh
LR::LRSplineSurface* ans = lrspline->copy(); LR::LRSplineSurface* ans = lrspline->copy();
ans->rebuildDimension(B.cols()); ans->rebuildDimension(points.rows());
// Swap around the control point values // Back to interleaved data
for (size_t j = 0; j < B.cols(); j++) { std::vector<double> interleave;
size_t i = 0; interleave.reserve(B.dim());
for (LR::Basisfunction* b : ans->getAllBasisfunctions()) for (size_t i = 0; i < nBasis; ++i)
b->cp()[j] = B(++i,j+1); for (size_t j = 0; j < points.rows(); j++) {
interleave.push_back(B(1+j*points.cols()+i));
} }
ans->setControlPoints(interleave);
return ans; return ans;
} }