diff --git a/CHANGELOG.md b/CHANGELOG.md
index 87bb694078c4ada30103e76ec4f04a3fd8a74126..8ebe450ad46632a25bf64de707f75679c9ffbb5c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,7 @@
# Master (will become release 2.8)
+- `UMFPackSolver` can now handle matrices and vectors with scalar entries.
+
- CholmodSolver: Added `errorCode_` member to report errors and warnings during the matrix decomposition
## Deprecations and removals
diff --git a/dune/solvers/common/defaultbitvector.hh b/dune/solvers/common/defaultbitvector.hh
index c8f301148d7ce8d134929eabccc74d1de80742f8..0b7d7cf3c26c86c4d616b9c559d4a7eac109e671 100644
--- a/dune/solvers/common/defaultbitvector.hh
+++ b/dune/solvers/common/defaultbitvector.hh
@@ -18,11 +18,18 @@ namespace Solvers {
namespace Imp {
+/** \brief Construct a bitvector type with a nesting that matches the one of a given vector
+ *
+ * The default implementation implements the case where the vector is a number type
+ * (to end recursion), and returns 'char' in such cases.
+ *
+ * \tparam Vector The vector types whose nesting shall be matched
+ */
template<class Vector>
struct DefaultBitVector
{
- static_assert(AlwaysFalse<Vector>::value, "No DefaultBitVector known for this type.");
- using type = void;
+ static_assert(IsNumber<Vector>::value, "No DefaultBitVector known for this type.");
+ using type = char;
};
template<class T, int i>
diff --git a/dune/solvers/norms/twonorm.hh b/dune/solvers/norms/twonorm.hh
index a62d6586d4bb1bf85da0c58e95500fa4111200be..ad06cde33637aed8cc7dbb7667a7b41db197bc9a 100644
--- a/dune/solvers/norms/twonorm.hh
+++ b/dune/solvers/norms/twonorm.hh
@@ -50,7 +50,7 @@ class TwoNorm : public Norm<V>
typename VectorType::block_type block = f1[i];
block -= f2[i];
- r += block.two_norm2();
+ r += Dune::Impl::asVector(block).two_norm2();
}
return std::sqrt(r);
diff --git a/dune/solvers/solvers/umfpacksolver.hh b/dune/solvers/solvers/umfpacksolver.hh
index 6071309190abf86846cadbed76f6013024fbf910..195d632cdca65ef41374d9e7d8ce018db4693d29 100644
--- a/dune/solvers/solvers/umfpacksolver.hh
+++ b/dune/solvers/solvers/umfpacksolver.hh
@@ -88,10 +88,16 @@ public:
for (size_t i=0; i<matrix_->N(); i++)
{
auto const &ignore = this->ignore()[i];
- if (ignore.none())
- nonIgnoreRows.insert(i);
- else if (not ignore.all())
- DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
+ if constexpr (std::is_convertible<decltype(ignore), bool>::value)
+ {
+ if (!ignore)
+ nonIgnoreRows.insert(i);
+ } else {
+ if (ignore.none())
+ nonIgnoreRows.insert(i);
+ else if (not ignore.all())
+ DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
+ }
}
// Construct the solver
@@ -111,15 +117,25 @@ public:
shortRowCount = 0;
for (size_t i=0; i<matrix_->N(); i++)
{
- if (this->ignore()[i].all())
- continue;
+ if constexpr (std::is_convertible<decltype(this->ignore()[i]), const bool>::value) {
+ if (this->ignore()[i])
+ continue;
+ } else {
+ if (this->ignore()[i].all())
+ continue;
+ }
auto cIt = (*matrix_)[i].begin();
auto cEndIt = (*matrix_)[i].end();
for (; cIt!=cEndIt; ++cIt)
- if (this->ignore()[cIt.index()].all())
- cIt->mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+ if constexpr (std::is_convertible<decltype(this->ignore()[cIt.index()]), const bool>::value) {
+ if (this->ignore()[cIt.index()])
+ Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+ } else {
+ if (this->ignore()[cIt.index()].all())
+ Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+ }
shortRowCount++;
}
diff --git a/dune/solvers/test/common.hh b/dune/solvers/test/common.hh
index 7015536debd45974e770bb404e2b8fed4f9c5275..62c48754495926ae80e5d78ca825a545336a8914 100644
--- a/dune/solvers/test/common.hh
+++ b/dune/solvers/test/common.hh
@@ -18,7 +18,7 @@
#include <dune/grid/utility/structuredgridfactory.hh>
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
#include <dune/grid/uggrid.hh>
#endif
@@ -26,6 +26,9 @@
#include <dune/alugrid/grid.hh>
#endif
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh>
@@ -64,13 +67,6 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
indices.exportIdx(matrix);
}
-template<class T, int n, int m, class F>
-void addToDiagonal(Dune::FieldMatrix<T, n, m>& M, const F& x)
-{
- for(int i=0; i<n; ++i)
- M[i][i] += x;
-}
-
template<class GridView, class Matrix>
void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
{
@@ -131,11 +127,11 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
double zij = (gradients[i] * gradients[j]) * z;
- addToDiagonal(matrix[iGlobal][jGlobal], zij);
- addToDiagonal(matrix[jGlobal][iGlobal], zij);
+ Dune::MatrixVector::addToDiagonal(matrix[iGlobal][jGlobal], zij);
+ Dune::MatrixVector::addToDiagonal(matrix[jGlobal][iGlobal], zij);
}
double zii = (gradients[i] * gradients[i]) * z;
- addToDiagonal(matrix[iGlobal][iGlobal], zii);
+ Dune::MatrixVector::addToDiagonal(matrix[iGlobal][iGlobal], zii);
}
}
}
@@ -314,15 +310,19 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
}
}
+/**
+ *
+ * \tparam blocksize If this is zero, then double will be used instead of FieldVector/FieldMatrix
+ */
template <size_t blocksize, class GridView, bool trivialDirichletOnly = true,
bool makePositiveDefinitive = not trivialDirichletOnly>
class SymmetricSampleProblem {
public:
- typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
- typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
+ using MatrixBlock = std::conditional_t<blocksize==0, double, Dune::FieldMatrix<double, blocksize, blocksize> >;
+ using VectorBlock = std::conditional_t<blocksize==0, double, Dune::FieldVector<double, blocksize> >;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector;
- typedef typename Dune::BitSetVector<blocksize> BitVector;
+ using BitVector = Dune::Solvers::DefaultBitVector_t<Vector>;
typedef EnergyNorm<Matrix, Vector> Norm;
SymmetricSampleProblem(GridView const& gridView) {
@@ -333,15 +333,21 @@ public:
energyNorm.setMatrix(&A);
ignore.resize(A.N());
- ignore.unsetAll();
+ if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
+ ignore.unsetAll();
+ else
+ std::fill(ignore.begin(), ignore.end(), false);
if (trivialDirichletOnly)
markBoundaryDOFs(gridView, ignore);
else {
// Mark the first component only
- Dune::BitSetVector<blocksize> boundaryNodes(A.N(), false);
+ BitVector boundaryNodes(A.N(), false);
markBoundaryDOFs(gridView, boundaryNodes);
for (size_t i=0; i < boundaryNodes.size(); ++i)
- ignore[i][0] = boundaryNodes[i][0];
+ if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
+ ignore[i][0] = boundaryNodes[i][0];
+ else
+ ignore[i] = boundaryNodes[i];
}
u.resize(A.N());
@@ -357,10 +363,17 @@ public:
// Set up a random 'solution'
u = 0;
for (size_t i = 0; i < u.size(); ++i)
- for (size_t j = 0; j < blocksize; ++j) {
+ if constexpr (blocksize==0)
+ {
if (makePositiveDefinitive)
- A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
- u_ex[i][j] = distribution(generator);
+ A[i][i] += 0.5*std::abs(A[0][0]);
+ u_ex[i] = distribution(generator);
+ } else {
+ for (size_t j = 0; j < blocksize; ++j) {
+ if (makePositiveDefinitive)
+ A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
+ u_ex[i][j] = distribution(generator);
+ }
}
// Construct right hand side corresponding to the 'solution'
@@ -368,9 +381,15 @@ public:
// set dirichlet values
for (size_t i = 0; i < u.size(); ++i)
- for (size_t j = 0; j < blocksize; ++j)
- if (ignore[i][j])
- u[i][j] = u_ex[i][j];
+ if constexpr (blocksize==0)
+ {
+ if (ignore[i])
+ u[i] = u_ex[i];
+ } else {
+ for (size_t j = 0; j < blocksize; ++j)
+ if (ignore[i][j])
+ u[i][j] = u_ex[i][j];
+ }
}
Matrix A;
diff --git a/dune/solvers/test/gssteptest.cc b/dune/solvers/test/gssteptest.cc
index 424ea189a240bf2e010a980a122aaa220fb285b4..9867e06525fe4c078e73d1e414fcaa7ce0423aba 100644
--- a/dune/solvers/test/gssteptest.cc
+++ b/dune/solvers/test/gssteptest.cc
@@ -250,7 +250,7 @@ bool checkWithSmallGrids(TestSuite& suite)
passed = passed and checkWithGrid<Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double,2> > >(suite, 3);
passed = passed and checkWithGrid<Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double,3> > >(suite, 1);
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
passed = passed and checkWithGrid<Dune::UGGrid<2> >(suite, 2);
passed = passed and checkWithGrid<Dune::UGGrid<3> >(suite, 1);
#endif
diff --git a/dune/solvers/test/umfpacksolvertest.cc b/dune/solvers/test/umfpacksolvertest.cc
index 545170c12951d75ecc147b95b3e04ea3027d8cc5..c396b26aa60dd68fc6ffeb903b364e387dd767e1 100644
--- a/dune/solvers/test/umfpacksolvertest.cc
+++ b/dune/solvers/test/umfpacksolvertest.cc
@@ -68,8 +68,10 @@ int main(int argc, char** argv)
bool passed = true;
+ UMFPackSolverTestSuite<0> testsuite0;
UMFPackSolverTestSuite<1> testsuite1;
UMFPackSolverTestSuite<2> testsuite2;
+ passed = checkWithStandardGrids(testsuite0);
passed = checkWithStandardGrids(testsuite1);
passed = checkWithStandardGrids(testsuite2);