localinterpolationconstraints.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_LOCALINTERPOLATIONCONSTRAINTS_HH
#define DUNE_FUFEM_CONSTRAINTS_LOCALINTERPOLATIONCONSTRAINTS_HH
 
#include <cmath>
#include <type_traits>
#include <vector>
 
#include <dune/common/rangeutilities.hh>
#include <dune/common/typetree/childaccess.hh>
#include <dune/common/typetree/traversal.hh>
 
#include <dune/functions/backends/istlvectorbackend.hh>
 
 
 
namespace Dune::Fufem {
 
 
 
template<class PrimaryLocalView, class SecondaryLocalView, class SecondaryToPrimary, class SecondaryDOFs, class Constraints>
void computeLocalInterpolationConstraints(const PrimaryLocalView& primaryLocalView, const SecondaryLocalView& secondaryLocalView, const SecondaryToPrimary& localSecondaryToPrimary, const SecondaryDOFs& secondaryDOFs, Constraints& constraints)
{
  using Coefficient = typename Constraints::Coefficient;
 
  auto&& isInterpolated = Dune::Functions::istlVectorBackend(constraints.isConstrained());
  auto&& interpolation = constraints.linearTerm();
 
  // Interpolation weights smaller than the tolerance will be ignored
  const auto tol = 1e-5;
 
  auto localIsConstrained = std::vector<bool>();
  auto interpolationValues = std::vector<Coefficient>();
 
  Dune::TypeTree::forEachLeafNode(secondaryLocalView.tree(), [&](auto&& secondaryNode, auto&& treePath) {
    auto&& primaryNode = Dune::TypeTree::child(primaryLocalView.tree(), treePath);
    if ((primaryNode.size() == 0) or (secondaryNode.size() == 0))
      return;
 
    const auto& primaryBasis = primaryNode.finiteElement().localBasis();
    using Range = typename std::decay_t<decltype(primaryBasis)>::Traits::RangeType;
 
    interpolationValues.resize(secondaryNode.finiteElement().size());
    localIsConstrained.clear();
    localIsConstrained.resize(secondaryNode.finiteElement().size(), false);
 
    // We first compute the primary basis values at all interpolation points.
    auto primaryBasisValues = std::vector<std::vector<Range>>(primaryBasis.size());
    auto k = 0;
    auto trackingBasisFunction = [&](const auto& xLocalSecondary) {
      primaryBasisValues.resize(k+1);
      auto xLocalPrimary = localSecondaryToPrimary(xLocalSecondary);
      primaryBasis.evaluateFunction(xLocalPrimary, primaryBasisValues[k]);
      return primaryBasisValues[k++][0];
    };
    secondaryNode.finiteElement().localInterpolation().interpolate(trackingBasisFunction, interpolationValues);
 
    // Node we interpolate any primary basis function on the current secondary element
    for(auto j : Dune::range(primaryBasis.size()))
    {
      auto j_global = primaryLocalView.index(primaryNode.localIndex(j));
 
      // primaryBasisFunction_j behaves like the j-th primary bases function if
      // interpolate() evaluates it at the same points in the same
      // order as above.
      auto k = 0;
      auto primaryBasisFunction_j = [&](const auto& x) {
        return primaryBasisValues[k++][j];
      };
      secondaryNode.finiteElement().localInterpolation().interpolate(primaryBasisFunction_j, interpolationValues);
 
      for(auto i : Dune::range(secondaryNode.finiteElement().size()))
      {
        // Only basis functions from the given set are considered
        // for being interpolated from the primary basis.
        if (not(secondaryDOFs.contains(secondaryNode.localIndex(i))))
          continue;
 
        auto i_global = secondaryLocalView.index(secondaryNode.localIndex(i));
 
        // We need to skip this if dof is already constrained.
        // Otherwise we may add interpolation weights from multiple elements
        // resulting in doubled interpolation weights after resolving the dependencies.
        if (isInterpolated[i_global])
          continue;
 
        if (i_global != j_global)
          if (std::abs(interpolationValues[i]) > tol)
          {
            interpolation[i_global][j_global] = interpolationValues[i];
            localIsConstrained[i] = true;
          }
      }
    }
    // Now mark all newly constrained nodes
    for(auto i : Dune::range(secondaryNode.finiteElement().size()))
      if (localIsConstrained[i])
        isInterpolated[secondaryLocalView.index(secondaryNode.localIndex(i))] = true;
  });
}
 
 
 
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_CONSTRAINTS_LOCALINTERPOLATIONCONSTRAINTS_HH