Skip to content
Snippets Groups Projects
Commit beecd24c authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Allow to generate test problems with scalar matrix entries

Previously, the class SymmetricSampleProblem would always produce
matrices and vectors with FieldMatrix/FieldVector entries of a given
user-controllable blocksize.  With this patch, the blocksize can
also be set to 0.  In that case, the resulting matrices and vectors
have 'double' entries.
parent 05d2f752
No related branches found
No related tags found
1 merge request!51UMFPackSolver with scalar matrix entries
...@@ -26,6 +26,9 @@ ...@@ -26,6 +26,9 @@
#include <dune/alugrid/grid.hh> #include <dune/alugrid/grid.hh>
#endif #endif
#include <dune/matrix-vector/addtodiagonal.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh> #include <dune/solvers/norms/twonorm.hh>
...@@ -64,13 +67,6 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix) ...@@ -64,13 +67,6 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
indices.exportIdx(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> template<class GridView, class Matrix>
void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix) void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
{ {
...@@ -131,11 +127,11 @@ 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); int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
double zij = (gradients[i] * gradients[j]) * z; double zij = (gradients[i] * gradients[j]) * z;
addToDiagonal(matrix[iGlobal][jGlobal], zij); Dune::MatrixVector::addToDiagonal(matrix[iGlobal][jGlobal], zij);
addToDiagonal(matrix[jGlobal][iGlobal], zij); Dune::MatrixVector::addToDiagonal(matrix[jGlobal][iGlobal], zij);
} }
double zii = (gradients[i] * gradients[i]) * z; 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) ...@@ -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, template <size_t blocksize, class GridView, bool trivialDirichletOnly = true,
bool makePositiveDefinitive = not trivialDirichletOnly> bool makePositiveDefinitive = not trivialDirichletOnly>
class SymmetricSampleProblem { class SymmetricSampleProblem {
public: public:
typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock; using MatrixBlock = std::conditional_t<blocksize==0, double, Dune::FieldMatrix<double, blocksize, blocksize> >;
typedef typename Dune::FieldVector<double, blocksize> VectorBlock; using VectorBlock = std::conditional_t<blocksize==0, double, Dune::FieldVector<double, blocksize> >;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix; typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector; typedef typename Dune::BlockVector<VectorBlock> Vector;
typedef typename Dune::BitSetVector<blocksize> BitVector; using BitVector = Dune::Solvers::DefaultBitVector_t<Vector>;
typedef EnergyNorm<Matrix, Vector> Norm; typedef EnergyNorm<Matrix, Vector> Norm;
SymmetricSampleProblem(GridView const& gridView) { SymmetricSampleProblem(GridView const& gridView) {
...@@ -333,15 +333,21 @@ public: ...@@ -333,15 +333,21 @@ public:
energyNorm.setMatrix(&A); energyNorm.setMatrix(&A);
ignore.resize(A.N()); ignore.resize(A.N());
if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
ignore.unsetAll(); ignore.unsetAll();
else
std::fill(ignore.begin(), ignore.end(), false);
if (trivialDirichletOnly) if (trivialDirichletOnly)
markBoundaryDOFs(gridView, ignore); markBoundaryDOFs(gridView, ignore);
else { else {
// Mark the first component only // Mark the first component only
Dune::BitSetVector<blocksize> boundaryNodes(A.N(), false); BitVector boundaryNodes(A.N(), false);
markBoundaryDOFs(gridView, boundaryNodes); markBoundaryDOFs(gridView, boundaryNodes);
for (size_t i=0; i < boundaryNodes.size(); ++i) for (size_t i=0; i < boundaryNodes.size(); ++i)
if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
ignore[i][0] = boundaryNodes[i][0]; ignore[i][0] = boundaryNodes[i][0];
else
ignore[i] = boundaryNodes[i];
} }
u.resize(A.N()); u.resize(A.N());
...@@ -357,21 +363,34 @@ public: ...@@ -357,21 +363,34 @@ public:
// Set up a random 'solution' // Set up a random 'solution'
u = 0; u = 0;
for (size_t i = 0; i < u.size(); ++i) for (size_t i = 0; i < u.size(); ++i)
if constexpr (blocksize==0)
{
if (makePositiveDefinitive)
A[i][i] += 0.5*std::abs(A[0][0]);
u_ex[i] = distribution(generator);
} else {
for (size_t j = 0; j < blocksize; ++j) { for (size_t j = 0; j < blocksize; ++j) {
if (makePositiveDefinitive) if (makePositiveDefinitive)
A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
u_ex[i][j] = distribution(generator); u_ex[i][j] = distribution(generator);
} }
}
// Construct right hand side corresponding to the 'solution' // Construct right hand side corresponding to the 'solution'
A.mv(u_ex, rhs); A.mv(u_ex, rhs);
// set dirichlet values // set dirichlet values
for (size_t i = 0; i < u.size(); ++i) for (size_t i = 0; i < u.size(); ++i)
if constexpr (blocksize==0)
{
if (ignore[i])
u[i] = u_ex[i];
} else {
for (size_t j = 0; j < blocksize; ++j) for (size_t j = 0; j < blocksize; ++j)
if (ignore[i][j]) if (ignore[i][j])
u[i][j] = u_ex[i][j]; u[i][j] = u_ex[i][j];
} }
}
Matrix A; Matrix A;
Vector u; Vector u;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment