symmetrize_constraints.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_LINEAR_SYMMETRIZE_CONSTRAINTS_HH
#define DUMUX_LINEAR_SYMMETRIZE_CONSTRAINTS_HH
 
#include <iostream>
 
#include <dune/common/hybridutilities.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/multitypeblockmatrix.hh>
 
namespace Dumux::Detail {
 
// loop over dense matrix block: end of recursion (actually do some work)
template<class K, int rows, int cols, class VectorRow, class VectorCol>
void symmetrizeConstraintsImpl(Dune::FieldMatrix<K, rows, cols>& A,
                               VectorRow& bRow, const VectorCol& bCol,
                               const VectorCol& constrainedRows,
                               bool isDiagonal)
{
    for (int j = 0; j < A.M(); ++j)
    {
        // all rows in this column can be eliminated by a Dirichlet row
        const auto& constrainedRowsJ = constrainedRows[j];
        if (constrainedRowsJ > 0.5)
        {
            for (int i = 0; i < A.N(); ++i)
            {
                auto& Aij = A[i][j];
                auto& bi = bRow[i];
                const auto& bj = bCol[j];
                if (!((i == j) && isDiagonal))
                {
                    bi -= Aij*bj;
                    Aij = 0.0;
                }
            }
        }
    }
}
 
// recursively loop over sub-blocks
template<class MBlock, class VectorRow, class VectorCol>
void symmetrizeConstraintsImpl(Dune::BCRSMatrix<MBlock>& A,
                               VectorRow& bRow, const VectorCol& bCol,
                               const VectorCol& constrainedRows,
                               bool isDiagonal)
{
    const auto rowEnd = A.end();
    for (auto row = A.begin(); row != rowEnd; ++row)
    {
        const auto colEnd = row->end();
        for (auto col = row->begin(); col != colEnd; ++col)
        {
            const auto i = row.index();
            const auto j = col.index();
 
            auto& Aij = A[i][j];
            auto& bi = bRow[i];
            const auto& bj = bCol[j];
            const auto& dj = constrainedRows[j];
 
            symmetrizeConstraintsImpl(Aij, bi, bj, dj, ((i == j) && isDiagonal));
        }
    }
}
 
// recursively loop over sub-blocks
template<class... MBlock, class VectorRow, class VectorCol>
void symmetrizeConstraintsImpl(Dune::MultiTypeBlockMatrix<MBlock...>& A,
                               VectorRow& bRow, const VectorCol& bCol,
                               const VectorCol& constrainedRows,
                               bool isDiagonal)
{
    using namespace Dune::Hybrid;
    forEach(std::make_index_sequence<Dune::MultiTypeBlockMatrix<MBlock...>::N()>{}, [&](const auto i)
    {
        forEach(std::make_index_sequence<Dune::MultiTypeBlockMatrix<MBlock...>::M()>{}, [&](const auto j)
        {
            auto& Aij = A[i][j];
            auto& bi = bRow[i];
            const auto& bj = bCol[j];
            const auto& dj = constrainedRows[j];
 
            symmetrizeConstraintsImpl(Aij, bi, bj, dj, ((i == j) && isDiagonal));
        });
    });
}
 
} // end namespace Dumux::Detail
 
namespace Dumux {
 
template<class MBlock, class Vector>
void symmetrizeConstraints(Dune::BCRSMatrix<MBlock>& A, Vector& b, const Vector& constrainedRows)
{
    Detail::symmetrizeConstraintsImpl(A, b, b, constrainedRows, true);
}
 
template<class... MBlock, class Vector>
void symmetrizeConstraints(Dune::MultiTypeBlockMatrix<MBlock...>& A, Vector& b, const Vector& constrainedRows)
{
    Detail::symmetrizeConstraintsImpl(A, b, b, constrainedRows, true);
}
 
} // end namespace Dumux
 
#endif