From b822c034ba95a75912e764bcb250581dcf9db50b Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Wed, 14 Nov 2018 10:22:39 +0100
Subject: [PATCH] [linear] Add functionality to omit values under threshold
Per default nothing is omitted. When a positive non-zero threshold is set,
all matrix entries below that threshold will be removed from the
pattern.
---
dumux/linear/matrixconverter.hh | 13 +++++++++++--
1 file changed, 11 insertions(+), 2 deletions(-)
diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh
index 39c2838608..aa15e391c3 100644
--- a/dumux/linear/matrixconverter.hh
+++ b/dumux/linear/matrixconverter.hh
@@ -24,6 +24,7 @@
#ifndef DUMUX_MATRIX_CONVERTER
#define DUMUX_MATRIX_CONVERTER
+#include <cmath>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/bvector.hh>
@@ -86,6 +87,9 @@ private:
// lambda function to fill the occupation pattern
auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
{
+ using std::abs;
+ static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
+
using BlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
const auto blockSizeI = BlockType::rows;
const auto blockSizeJ = BlockType::cols;
@@ -93,7 +97,8 @@ private:
for(auto col = row->begin(); col != row->end(); ++col)
for(std::size_t i = 0; i < blockSizeI; ++i)
for(std::size_t j = 0; j < blockSizeJ; ++j)
- occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j);
+ if(abs(subMatrix[row.index()][col.index()][i][j]) > eps)
+ occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j);
};
// fill the pattern
@@ -132,6 +137,9 @@ private:
// lambda function to copy the values
auto copyValues = [&M](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
{
+ using std::abs;
+ static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
+
using BlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
const auto blockSizeI = BlockType::rows;
const auto blockSizeJ = BlockType::cols;
@@ -139,7 +147,8 @@ private:
for (auto col = row->begin(); col != row->end(); ++col)
for (std::size_t i = 0; i < blockSizeI; ++i)
for (std::size_t j = 0; j < blockSizeJ; ++j)
- M[startRow + row.index()*blockSizeI + i][startCol + col.index()*blockSizeJ + j] = subMatrix[row.index()][col.index()][i][j];
+ if(abs(subMatrix[row.index()][col.index()][i][j]) > eps)
+ M[startRow + row.index()*blockSizeI + i][startCol + col.index()*blockSizeJ + j] = subMatrix[row.index()][col.index()][i][j];
};
--
GitLab