continuityconstraints.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_CONTINUITYCONSTRAINTS_HH
#define DUNE_FUFEM_CONSTRAINTS_CONTINUITYCONSTRAINTS_HH
 
#include <type_traits>
#include <variant>
 
#include <dune/common/rangeutilities.hh>
 
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/subentitydofs.hh>
 
#include <dune/fufem/constraints/affineconstraints.hh>
#include <dune/fufem/constraints/localinterpolationconstraints.hh>
 
 
 
namespace Dune::Fufem {
 
 
 
template<class BV, class V, class MI, class C, class Basis>
void computeContinuityConstraints(AffineConstraints<BV, V, MI, C>& constraints, const Basis& basis)
{
  using Constraints = AffineConstraints<BV, V, MI, C>;
 
  static_assert(std::is_same_v<typename Basis::MultiIndex, typename Constraints::MultiIndex>);
 
  auto&& isInterpolated = Dune::Functions::istlVectorBackend(constraints.isConstrained());
 
  auto localIndices = [](const auto& tree) {
    if (tree.size()==0)
      return Dune::range(std::size_t(0), std::size_t(0));
    return Dune::range(tree.localIndex(0), tree.localIndex(0)+tree.size());
  };
 
  auto localView = basis.localView();
  auto neighborLocalView = basis.localView();
  auto intersectionDOFs = Dune::Functions::subEntityDOFs(basis);
 
  for(const auto& element : elements(basis.gridView()))
  {
    localView.bind(element);
    auto level = element.level();
 
    bool allDofsConstrained = true;
    for(auto i : localIndices(localView.tree()))
      allDofsConstrained = allDofsConstrained and isInterpolated[localView.index(i)];
    if (allDofsConstrained)
        continue;
 
    for(const auto& intersection : intersections(basis.gridView(), element))
    {
      if (not(intersection.neighbor()))
        continue;
 
      const auto& neighborElement = intersection.outside();
 
      if (neighborElement.level()>=level)
        continue;
 
      intersectionDOFs.bind(localView, intersection);
 
      neighborLocalView.bind(neighborElement);
 
      auto secondaryToPrimary = [&](const auto& x) {
        auto globalPosition = element.geometry().global(x);
        return neighborElement.geometry().local(globalPosition);
      };
 
      computeLocalInterpolationConstraints(neighborLocalView, localView, secondaryToPrimary, intersectionDOFs, constraints);
    }
  }
  constraints.resolveDependencies();
}
 
 
 
template<class BitVector=std::monostate, class Vector=std::monostate, class Basis>
auto makeContinuityConstraints(const Basis& basis)
{
  auto constraints = makeAffineConstraints<BitVector, Vector>(basis);
  computeContinuityConstraints(constraints, basis);
  return constraints;
}
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_CONSTRAINTS_CONTINUITYCONSTRAINTS_HH