affineconstraints.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_CONSTRAINTS_AFFINECONSTRAINTS_HH
#define DUNE_FUFEM_CONSTRAINTS_AFFINECONSTRAINTS_HH
 
#include <iostream>
#include <map>
#include <unordered_map>
#include <utility>
#include <type_traits>
 
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/backends/containerfactory.hh>
#include <dune/functions/backends/istlvectorfactory.hh>
 
#include <dune/assembler/backends/istlmatrixbackend.hh>
 
#include <dune/fufem/backends/traversal.hh>
 
namespace Dune::Fufem {
 
namespace Impl {
 
  // This class stores a range of indices and wraps it in
  // a LocalView-like interface. This allows to pass a
  // generic set of indices to the scatter() methods
  // of the MatrixBackend implementations in dune-assembler.
  template<class Index>
  class FakeLocalView
  {
    class FakeTree
    {
      std::size_t size_;
 
    public:
 
      FakeTree(std::size_t size)
        : size_(size)
      {}
 
      constexpr std::size_t size() const
      {
        return size_;
      }
 
      constexpr std::size_t localIndex(std::size_t nodeLocalIndex) const
      {
        return nodeLocalIndex;
      }
    };
 
    using Indices = std::vector<Index>;
 
    Indices indices_;
 
  public:
 
    constexpr FakeTree tree() const
    {
      return FakeTree{size()};
    }
 
    constexpr std::size_t size() const
    {
      return indices_.size();
    }
 
    constexpr const Index& index(std::size_t localIndex) const
    {
      return indices_[localIndex];
    }
 
    constexpr const auto& indices() const
    {
      return indices_;
    }
 
    constexpr auto& indices()
    {
      return indices_;
    }
 
    template<class Range>
    constexpr void fillFromRange(Range&& range)
    {
      indices_.resize(range.size());
      std::copy(range.begin(), range.end(), indices_.begin());
    }
 
    constexpr void clear()
    {
      indices_.clear();
    }
 
  };
 
} // end namespace Impl
 
 
 
template<class BV, class V, class MI, class C = double>
class AffineConstraints
{
 
  template<class OtherMultiIndex>
  static decltype(auto) toMultiIndex(const OtherMultiIndex& omi)
  {
    if constexpr (std::is_same_v<OtherMultiIndex, MultiIndex>)
      return omi;
    else if constexpr (std::is_convertible_v<OtherMultiIndex, MultiIndex>)
      return MultiIndex(omi);
    else if constexpr (requires { MultiIndex{{omi}}; })
      return MultiIndex{{omi}};
    else if constexpr (requires(MultiIndex mi){ mi.push_back(0); omi.size(); })
    {
      MultiIndex mi;
      for(auto i : omi)
        mi.push_back(i);
      return mi;
    }
    else
      static_assert(Dune::AlwaysFalse<OtherMultiIndex>::value, "Value cannot be converted to MultiIndex type");
  }
 
  // The following two helpers exist, because the vector backends support non-matching
  // container types, where a `FieldVector<T,1>` is used instead of a plain `T`. Then
  // traversing the container will produce a trailing multi-index digit 0 which does
  // not show up in the basis. Hence, a search in the unordered_map fails if the
  // multi-index has been inserted as key without the trailing 0. To avoid this pitfall
  // the following helpers ignore a trailing multi-index digit 0 when computing hashes
  // and during equality comparison. Notice that this does not produce false positives
  // because a basis cannot contain the multi-indices `a=(a_0,...,a_n,0)` and `b=(a_0,...a_n)`
  // at the same time, since the latter would correspond to an inner node in the index
  // tree if the former exists.
  //
  // These custom functors can be removed once we require that the container nesting
  // structure must exactly match the basis indices and check for this during index access.
  // Until then, using these functors avoids that some constraints are silently
  // overlooked, because they are not found by `std::unordered_map::find` due
  // to the mismatch.
  struct MultiIndexHash
  {
    std::size_t operator()(const MI& mi) const
    {
      std::size_t seed = 0;
      for (auto i : Dune::range(mi.size()-1))
        Dune::hash_combine(seed, mi[i]);
      if (mi.back() != 0)
        Dune::hash_combine(seed, mi.back());
      return seed;
    }
  };
 
  struct MultiIndexEquals
  {
    bool operator()(const MI& mi1, const MI& mi2) const
    {
      auto compareDigits = [&](auto size) {
        for (auto i : Dune::range(size))
          if (mi1[i] != mi2[i])
            return false;
        return true;
      };
      if (mi1.size() == mi2.size())
        return compareDigits(mi1.size());
      else if (mi1.size()+1 == mi2.size())
      {
        if (mi2.back() == 0)
          return compareDigits(mi1.size());
      }
      else if (mi2.size()+1 == mi1.size())
      {
        if (mi1.back() == 0)
          return compareDigits(mi2.size());
      }
      return false;
    }
  };
 
public:
 
  using BitVector = BV;
  using Vector = V;
  using MultiIndex = MI;
  using Coefficient = C;
 
  using LinearTermRow = std::map<MI, Coefficient>;
  using LinearTerm = std::unordered_map<MultiIndex, LinearTermRow, MultiIndexHash, MultiIndexEquals>;
  using ConstantTerm = Vector;
 
 
 
  AffineConstraints() = default;
  AffineConstraints(const AffineConstraints&) = default;
  AffineConstraints(AffineConstraints&&) = default;
 
  AffineConstraints(BitVector&& isConstrained, Vector&& constantTerm)
    : isConstrained_(std::move(isConstrained))
    , constantTerm_(std::move(constantTerm))
  {}
 
  AffineConstraints(const BitVector& isConstrained, const Vector& constantTerm)
    : isConstrained_(isConstrained)
    , constantTerm_(constantTerm)
  {}
 
  template<class MultiIndexType>
  bool isConstrained(const MultiIndexType& i) const
  {
    return Dune::Functions::istlVectorBackend(isConstrained_)[toMultiIndex(i)];
  }
 
  const BitVector& isConstrained() const
  {
    return isConstrained_;
  }
 
  BitVector& isConstrained()
  {
    return isConstrained_;
  }
 
  const LinearTerm& linearTerm() const
  {
    return interpolation_;
  }
 
  LinearTerm& linearTerm()
  {
    return interpolation_;
  }
 
  template<class MultiIndexType>
  const LinearTermRow& linearTerm(const MultiIndexType& i) const
  {
    auto it = interpolation_.find(toMultiIndex(i));
    if (it != interpolation_.end())
      return it->second;
    else
      return emptyInterpolationRow_;
  }
 
  const ConstantTerm& constantTerm() const
  {
    return constantTerm_;
  }
 
  ConstantTerm& constantTerm()
  {
    return constantTerm_;
  }
 
  template<class MultiIndexType>
  const auto& constantTerm(const MultiIndexType& i) const
  {
    return Dune::Functions::istlVectorBackend(constantTerm_)[toMultiIndex(i)];
  }
 
  template<class MultiIndexType>
  auto& constantTerm(const MultiIndexType& i)
  {
    return Dune::Functions::istlVectorBackend(constantTerm_)[i];
  }
 
  template<class Vector>
  void interpolate(Vector& v) const
  {
    auto vectorBackend = Dune::Functions::istlVectorBackend(v);
    Dune::Fufem::recursiveForEachVectorEntry(v, [&](auto&& vector_i, const auto& i)
    {
      if (isConstrained(i))
      {
        vector_i = constantTerm(i);
        for(const auto& [j, L_ij] : linearTerm(i))
          vector_i += L_ij * vectorBackend[j];
      }
    });
  }
 
  template<class Vector>
  void interpolateHomogeneously(Vector& v) const
  {
    auto vectorBackend = Dune::Functions::istlVectorBackend(v);
    Dune::Fufem::recursiveForEachVectorEntry(v, [&](auto&& vector_i, const auto& i)
    {
      if (isConstrained(i))
      {
        vector_i = 0;
        for(const auto& [j, L_ij] : linearTerm(i))
          vector_i += L_ij * vectorBackend[j];
      }
    });
  }
 
  template<class MatrixPatternBackend, class Basis>
  void extendMatrixPattern(MatrixPatternBackend& patternBackend, const Basis& basis) const
  {
    // Check if the MatrixPatternBackend supports the dune-assembler interface.
    if constexpr (requires() { typename MatrixPatternBackend::LocalPattern; })
    {
      // This functions extracts the global indices from a linear combination
      auto globalIndices = [](const auto& linearTerm) {
        return Dune::transformedRangeView(
          linearTerm,
          [](const auto& coeff) -> decltype(auto) { return coeff.first; }
        );
      };
 
      // We need to pass LocalViews and a LocalPattern to the scatter
      // methods of dune-assembler MatrixBackends.
      auto localPattern = typename MatrixPatternBackend::LocalPattern();
      auto newRowIndices = Impl::FakeLocalView<typename Basis::MultiIndex>();
      auto newColIndices = Impl::FakeLocalView<typename Basis::MultiIndex>();
      auto accumulatedRow = typename Basis::MultiIndex();
      auto accumulatedColIndices = std::vector<typename Basis::MultiIndex>();
      patternBackend.forEachEntry([&](const auto& row, const auto& col)
      {
        if (isConstrained(row))
        {
          newRowIndices.fillFromRange(globalIndices(linearTerm(toMultiIndex(row))));
          if (isConstrained(col))
            newColIndices.fillFromRange(globalIndices(linearTerm(toMultiIndex(col))));
          else
            newColIndices.fillFromRange(std::array{toMultiIndex(col)});
          localPattern.resize(newRowIndices.size(), newColIndices.size());
          localPattern.addAll();
          patternBackend.scatter(newRowIndices, newColIndices, localPattern);
        }
        else
        {
          if (isConstrained(col))
          {
            if (accumulatedRow != toMultiIndex(row))
            {
              if (accumulatedColIndices.size()>0)
              {
                newRowIndices.fillFromRange(std::array{accumulatedRow});
                std::swap(newColIndices.indices(), accumulatedColIndices);
                localPattern.resize(newRowIndices.size(), newColIndices.size());
                localPattern.addAll();
                patternBackend.scatter(newRowIndices, newColIndices, localPattern);
              }
              accumulatedColIndices.clear();
              accumulatedRow = toMultiIndex(row);
            }
            for(const auto& newColIndex : globalIndices(linearTerm(toMultiIndex(col))))
              accumulatedColIndices.push_back(newColIndex);
          }
        }
      });
      if (accumulatedColIndices.size()>0)
      {
        newRowIndices.fillFromRange(std::array{accumulatedRow});
        std::swap(newColIndices.indices(), accumulatedColIndices);
        localPattern.resize(newRowIndices.size(), newColIndices.size());
        localPattern.addAll();
        patternBackend.scatter(newRowIndices, newColIndices, localPattern);
      }
    }
    else
    {
      auto localView = basis.localView();
      for(const auto& element : elements(basis.gridView()))
      {
        localView.bind(element);
        for(auto i_local : Dune::range(localView.tree().size()))
        {
          auto i = localView.index(localView.tree().localIndex(i_local));
          for(auto j_local : Dune::range(localView.tree().size()))
          {
            auto j = localView.index(localView.tree().localIndex(j_local));
            if (isConstrained(i))
            {
              if (isConstrained(j))
                for(const auto& [k, L_ik] : linearTerm(i))
                  for(const auto& [l, L_jl] : linearTerm(j))
                    patternBackend.insertEntry(k, l);
              else
                for(const auto& [k, L_ik] : linearTerm(i))
                  patternBackend.insertEntry(k, j);
            }
            else if (isConstrained(j))
            {
              for(const auto& [l, L_jl] : linearTerm(j))
                patternBackend.insertEntry(i, l);
            }
          }
        }
      }
    }
  }
 
  template<class Matrix>
  void constrainMatrix(Matrix& A) const
  {
    auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(A);
    Dune::Fufem::recursiveForEachMatrixEntry(A, [&](auto&& matrix_ij, const auto& i, const auto& j)
    {
      auto equals = [](const auto& i, const auto& j) {
        if constexpr (std::is_same_v<decltype(i), decltype(j)>)
          return i==j;
        else
          return false;
      };
      if (isConstrained(i))
      {
        if (isConstrained(j))
          for(const auto& [k, L_ik] : linearTerm(i))
            for(const auto& [l, L_jl] : linearTerm(j))
              matrixBackend(k, l) += L_ik*L_jl * matrix_ij;
        else
          for(const auto& [k, L_ik] : linearTerm(i))
            matrixBackend(k, j) += L_ik * matrix_ij;
        matrix_ij = equals(i, j) ? 1 : 0;
      }
      else if (isConstrained(j))
      {
        for(const auto& [l, L_jl] : linearTerm(j))
          matrixBackend(i, l) += L_jl * matrix_ij;
        matrix_ij = equals(i, j) ? 1 : 0;
      }
    });
  }
 
  template<class Vector>
  void constrainVector(Vector& b) const
  {
    auto vectorBackend = Dune::Functions::istlVectorBackend(b);
    Dune::Fufem::recursiveForEachVectorEntry(b, [&](auto&& vector_i, const auto& i)
    {
      if (isConstrained(i))
      {
        for(const auto& [k, L_ik] : linearTerm(i))
          vectorBackend[k] += L_ik * vector_i;
        vector_i = constantTerm(i);
      }
    });
  }
 
  template<class Vector>
  void constrainVectorHomogeneously(Vector& b) const
  {
    auto vectorBackend = Dune::Functions::istlVectorBackend(b);
    Dune::Fufem::recursiveForEachVectorEntry(b, [&](auto&& vector_i, const auto& i)
    {
      if (isConstrained(i))
      {
        for(const auto& [k, L_ik] : linearTerm(i))
          vectorBackend[k] += L_ik * vector_i;
        vector_i = 0;
      }
    });
  }
 
  template<class Matrix, class Vector>
  void constrainLinearSystem(Matrix& A, Vector& b) const
  {
    A.mmv(constantTerm_, b);
    constrainMatrix(A);
    constrainVector(b);
  }
 
  template<class Matrix, class Vector>
  void constrainLinearSystemHomogeneously(Matrix& A, Vector& b) const
  {
    constrainMatrix(A);
    constrainVectorHomogeneously(b);
  }
 
  void report() const
  {
    Dune::Fufem::recursiveForEachVectorEntry(isConstrained_, [&](auto&& isConstrained_i, const auto& i)
    {
      if (isConstrained(i))
      {
        std::cout << "x_" << i << " = ";
        for(const auto& [j, L_ij] : linearTerm(i))
          std::cout << L_ij << " * x_" << j << " + " ;
        std::cout << constantTerm(i) << std::endl;
      }
    });
    std::cout << "interpolated dofs : " << interpolation_.size() << std::endl;
  }
 
  bool check() const
  {
    bool pass = true;
    for(const auto& [i, L_i] : interpolation_)
    {
      for(const auto& [j, L_ij] : L_i)
      {
        if (isConstrained(j))
        {
          std::cout << "Error: interpolation for constrained DOF " << i << " depends on constrained DOF " << j << std::endl;
          pass = false;
        }
      }
    }
    return pass;
  }
 
  void resolveDependencies()
  {
    auto dependenciesResolved = isConstrained_;
    auto dependenciesResolvedBackend = Dune::Functions::istlVectorBackend(dependenciesResolved);
    dependenciesResolvedBackend = false;
    for(auto& [i, L_i] : interpolation_)
      if (not(dependenciesResolvedBackend[i]))
        resolveDOF(dependenciesResolved, i, L_i);
  }
 
private:
 
  void resolveDOF(BitVector& dependenciesResolved, const MultiIndex& i, LinearTermRow& L_i)
  {
    auto dependenciesResolvedBackend = Dune::Functions::istlVectorBackend(dependenciesResolved);
    auto c = Dune::Functions::istlVectorBackend(constantTerm_);
    auto it = L_i.begin();
    while (it != L_i.end())
    {
      auto& [j, L_ij] = *it;
      // If interpolation weight is itself constraint,
      // first resolve it recursively (if needed)
      // and then insert all indirect dependencies directly
      if (isConstrained(j))
      {
        c[i] += L_ij * c[j];
        auto L_j_it = interpolation_.find(j);
        if (L_j_it != interpolation_.end())
        {
          auto& L_j = L_j_it->second;
          if (not(dependenciesResolvedBackend[j]))
            resolveDOF(dependenciesResolved, j, L_j);
          for (const auto& [k, L_jk] : L_j)
            L_i[k] += L_ij * L_jk;
        }
        else
          dependenciesResolvedBackend[j] = true;
        it = L_i.erase(it);
      }
      else
        ++it;
    }
    dependenciesResolvedBackend[i] = true;
  }
 
  BitVector isConstrained_;
  LinearTerm interpolation_;
  LinearTermRow emptyInterpolationRow_;
  ConstantTerm constantTerm_;
};
 
 
 
template<class BitVector=std::monostate, class Vector=std::monostate, class Basis>
auto makeAffineConstraints(const Basis& basis)
{
  using MultiIndex = typename Basis::MultiIndex;
  using Coefficient = double;
 
  auto generateBitVector = [&] {
    if constexpr (std::is_same_v<BitVector, std::monostate>)
      return Dune::Functions::makeContainer<char>(basis.rootBasis().preBasis().containerDescriptor(), false);
    else
    {
      auto v = BitVector();
      Dune::Functions::istlVectorBackend(v).resize(basis);
      Dune::Functions::istlVectorBackend(v) = false;
      return v;
    }
  };
  auto generateVector = [&] {
    if constexpr (std::is_same_v<Vector, std::monostate>)
      return Dune::Functions::makeISTLVector<Coefficient>(basis.rootBasis().preBasis().containerDescriptor());
    else
    {
      auto v = Vector();
      Dune::Functions::istlVectorBackend(v).resize(basis);
      Dune::Functions::istlVectorBackend(v) = false;
      return v;
    }
  };
  using GeneratedBitVector = decltype(generateBitVector());
  using GeneratedVector = decltype(generateVector());
  static_assert(not std::is_same_v<GeneratedBitVector, std::monostate>, "Deducing a BitVector-type for non-flat indices in makeAffineConstraints requires dune-functions>=2.11.");
  static_assert(not std::is_same_v<GeneratedVector, std::monostate>, "Deducing a Vector-type for non-flat indices in makeAffineConstraints requires dune-functions>=2.11.");
  return AffineConstraints<GeneratedBitVector, GeneratedVector, MultiIndex, Coefficient>(generateBitVector(), generateVector());
}
 
 
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_CONSTRAINTS_AFFINECONSTRAINTS_HH