GCC Code Coverage Report

Line	Branch	Exec	Source
1			// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2			// vi: set et ts=4 sw=4 sts=4:
3			//
4			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5			// SPDX-License-Identifier: GPL-3.0-or-later
6			//
7			/*!
8			* \file
9			* \ingroup Linear
10			* \brief Helper function to symmetrize row-constraints in constrained linear systems
11			*/
12			#ifndef DUMUX_LINEAR_SYMMETRIZE_CONSTRAINTS_HH
13			#define DUMUX_LINEAR_SYMMETRIZE_CONSTRAINTS_HH
14
15			#include <iostream>
16
17			#include <dune/common/hybridutilities.hh>
18			#include <dune/common/fmatrix.hh>
19			#include <dune/istl/bcrsmatrix.hh>
20			#include <dune/istl/multitypeblockmatrix.hh>
21
22			namespace Dumux::Detail {
23
24			// loop over dense matrix block: end of recursion (actually do some work)
25			template<class K, int rows, int cols, class VectorRow, class VectorCol>
26			void symmetrizeConstraintsImpl(Dune::FieldMatrix<K, rows, cols>& A,
27			VectorRow& bRow, const VectorCol& bCol,
28			const VectorCol& constrainedRows,
29			bool isDiagonal)
30			{
31		✗	for (int j = 0; j < A.M(); ++j)
32			{
33			// all rows in this column can be eliminated by a Dirichlet row
34		✗	const auto& constrainedRowsJ = constrainedRows[j];
35		✗	if (constrainedRowsJ > 0.5)
36			{
37		✗	for (int i = 0; i < A.N(); ++i)
38			{
39		✗	auto& Aij = A[i][j];
40		✗	auto& bi = bRow[i];
41		✗	const auto& bj = bCol[j];
42		✗	if (!((i == j) && isDiagonal))
43			{
44		✗	bi -= Aij*bj;
45		✗	Aij = 0.0;
46			}
47			}
48			}
49			}
50			}
51
52			// recursively loop over sub-blocks
53			template<class MBlock, class VectorRow, class VectorCol>
54		✗	void symmetrizeConstraintsImpl(Dune::BCRSMatrix<MBlock>& A,
55			VectorRow& bRow, const VectorCol& bCol,
56			const VectorCol& constrainedRows,
57			bool isDiagonal)
58			{
59		✗	const auto rowEnd = A.end();
60		✗	for (auto row = A.begin(); row != rowEnd; ++row)
61			{
62		✗	const auto colEnd = row->end();
63		✗	for (auto col = row->begin(); col != colEnd; ++col)
64			{
65		✗	const auto i = row.index();
66		✗	const auto j = col.index();
67
68		✗	auto& Aij = A[i][j];
69		✗	auto& bi = bRow[i];
70		✗	const auto& bj = bCol[j];
71		✗	const auto& dj = constrainedRows[j];
72
73		✗	symmetrizeConstraintsImpl(Aij, bi, bj, dj, ((i == j) && isDiagonal));
74			}
75			}
76		✗	}
77
78			// recursively loop over sub-blocks
79			template<class... MBlock, class VectorRow, class VectorCol>
80		✗	void symmetrizeConstraintsImpl(Dune::MultiTypeBlockMatrix<MBlock...>& A,
81			VectorRow& bRow, const VectorCol& bCol,
82			const VectorCol& constrainedRows,
83			bool isDiagonal)
84			{
85			using namespace Dune::Hybrid;
86		✗	forEach(std::make_index_sequence<Dune::MultiTypeBlockMatrix<MBlock...>::N()>{}, [&](const auto i)
87			{
88		✗	forEach(std::make_index_sequence<Dune::MultiTypeBlockMatrix<MBlock...>::M()>{}, [&](const auto j)
89			{
90		✗	auto& Aij = A[i][j];
91		✗	auto& bi = bRow[i];
92		✗	const auto& bj = bCol[j];
93		✗	const auto& dj = constrainedRows[j];
94
95		✗	symmetrizeConstraintsImpl(Aij, bi, bj, dj, ((i == j) && isDiagonal));
96			});
97			});
98		✗	}
99
100			} // end namespace Dumux::Detail
101
102			namespace Dumux {
103
104			/*!
105			* \ingroup Linear
106			* \brief Symmetrize the constrained system Ax=b
107			* \note The function addresses the case where some rows in the linear system have been replaced
108			* to fix a degree of freedom to a fixed value. A typical use case is strongly enforced Dirichlet constraints.
109			* To symmetrize the constraints in operator A, this function eliminates the corresponding column entries by
110			* bringing them to the right-hand side. Particularly, if matrix operator A was symmetric before incorporating
111			* constraints, calling this function restores the symmetry of operator A.
112			* \note the constrainedRows vector is the same shape and type as b and contains 1 for each row
113			* that corresponds to constraints (e.g. strongly enforced Dirichlet dofs) and 0 elsewhere
114			* \note this is a specialization for Dune::BCRSMatrix
115			*/
116			template<class MBlock, class Vector>
117			void symmetrizeConstraints(Dune::BCRSMatrix<MBlock>& A, Vector& b, const Vector& constrainedRows)
118			{
119			Detail::symmetrizeConstraintsImpl(A, b, b, constrainedRows, true);
120			}
121
122			/*!
123			* \ingroup Linear
124			* \brief Symmetrize the constrained system Ax=b
125			* \note The function addresses the case where some rows in the linear system have been replaced
126			* to fix a degree of freedom to a fixed value. A typical use case is strongly enforced Dirichlet constraints.
127			* To symmetrize the constraints in operator A, this function eliminates the corresponding column entries by
128			* bringing them to the right-hand side. Particularly, if matrix operator A was symmetric before incorporating
129			* constraints, calling this function restores the symmetry of operator A.
130			* \note the constrainedRows vector is the same shape and type as b and contains 1 for each row
131			* that corresponds to constraints (e.g. strongly enforced Dirichlet dofs) and 0 elsewhere
132			* \note this is a specialization for Dune::MultiTypeBlockMatrix
133			*/
134			template<class... MBlock, class Vector>
135			void symmetrizeConstraints(Dune::MultiTypeBlockMatrix<MBlock...>& A, Vector& b, const Vector& constrainedRows)
136			{
137		✗	Detail::symmetrizeConstraintsImpl(A, b, b, constrainedRows, true);
138			}
139
140			} // end namespace Dumux
141
142			#endif
143