Commit 363c8174 authored by Lasse Hinrichsen's avatar Lasse Hinrichsen
Browse files

Generalize sumfunctional's linearization

parent b2f1c660
......@@ -2,13 +2,48 @@
#define DUNE_TNNMG_SUMFUNCTIONALCONSTRAINEDLINEARIZATION_HH
#include <tuple>
#include <type_traits>
#include <dune/common/concept.hh>
#include <dune/common/hybridutilities.hh>
// The following include is only to have the truncateMatrix/-Vector functions.
// These should be in their own file at some point.
#include <dune/tnnmg/functionals/semilinearfunctionallinearization.hh>
#include <dune/solvers/common/resize.hh>
namespace Dune::TNNMG
{
namespace Helper {
struct HasBitwiseOrAssignment
{
template<typename C, typename D>
auto require(C&& c, D&& d) -> decltype(c |= d);
};
/** \brief Compute |= on the bits of a possible nested bitvector.
*
* TODO: Reference corresponding MR.
*/
template<typename BV0, typename BV1>
void
bitwiseOrAssignment(BV0&& lhs, BV1&& rhs)
{
namespace H = Dune::Hybrid;
if constexpr (models < HasBitwiseOrAssignment,
decltype(lhs),
decltype(rhs)>()) {
lhs |= rhs;
} else {
H::forEach(H::integralRange(H::size(rhs)), [&](auto i) {
bitwiseOrAssignment(H::elementAt(lhs, i), H::elementAt(rhs, i));
});
}
}
}
/**
* \brief A constrained linearization for a sum of functionals
*
......@@ -37,6 +72,14 @@ public:
ignore_(ignore)
{}
/** \brief Constructor based on already constructed linearizations */
template<typename Linearizations>
SumFunctionalConstrainedLinearization(const Linearizations& linearizations,
const BitVector& ignore)
: addends_(linearizations)
, ignore_(ignore)
{}
/** \brief Pre-compute derivative information at the configuration x
*/
void bind(const Vector& x)
......@@ -67,37 +110,27 @@ public:
truncationFlags_ = ignore_;
// Also truncate a degree of freedom if one of the addends wants it truncated
for (std::size_t i=0; i<truncationFlags_.size(); i++)
Hybrid::forEach(addends_, [&](auto&& addend) {
truncationFlags_[i] |= addend.truncated()[i];
});
Hybrid::forEach(addends_, [&](auto&& addend) {
Helper::bitwiseOrAssignment(truncationFlags_, addend.truncated());
});
/////////////////////////////////////////////////////
// Do the actual truncation
/////////////////////////////////////////////////////
for (std::size_t i=0; i<x.size(); ++i)
{
// Truncate the gradient
for (std::size_t j=0; j<x[i].size(); ++j)
if (truncationFlags_[i][j])
negativeGradient_[i][j] = 0;
// Truncate the Hesse matrix
for (auto&& [entry, l] : sparseRange(hessian_[i]))
for (std::size_t j=0; j<x[i].size(); ++j)
for (std::size_t k=0; k<x[i].size(); ++k)
if (truncationFlags_[i][j] || truncationFlags_[l][k])
{
entry[j][k] = 0;
#ifndef GAUSSSEIDEL
if (j==k && i == l)
entry[j][k] = 1;
#endif
}
}
}
Helper::truncateVector(negativeGradient_, truncationFlags_);
Helper::truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
// TODO: Re-introduce this part for those poor people who need 1s on the
// truncated diagonal. However, this should probably be better documented
// than defining the GAUSSSEIDEL macro.
/*
#ifndef GAUSSSEIDEL
if (j==k && i == l)
entry[j][k] = 1;
#endif
*/
}
/** \brief Fill the truncated degrees of freedom in a vector with zeros
*
......@@ -107,9 +140,8 @@ public:
*/
void extendCorrection(const ConstrainedVector& cv, Vector& v) const
{
for (std::size_t i=0; i<v.size(); ++i)
for (std::size_t j=0; j<v[i].size(); ++j)
v[i][j] = (truncationFlags_[i][j]) ? 0 : cv[i][j];
v = cv;
Helper::truncateVector(v, truncationFlags_);
}
/** \brief Access to which degrees of freedom have been truncated */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment