periodicconstraints.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_PERIODICCONSTRAINTS_HH
#define DUNE_FUFEM_CONSTRAINTS_PERIODICCONSTRAINTS_HH
 
#include <type_traits>
#include <variant>
#include <vector>
#include <tuple>
 
#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 PrimaryBasis, class PrimaryPatch, class SecondaryBasis, class SecondaryPatch, class SecondaryToPrimary>
void computePeriodicConstraints(
  AffineConstraints<BV, V, MI, C>& constraints,
  const PrimaryBasis& primaryBasis,
  const PrimaryPatch& primaryPatch,
  const SecondaryBasis& secondaryBasis,
  const SecondaryPatch& secondaryPatch,
  const SecondaryToPrimary& secondaryToPrimary)
{
  using Constraints = AffineConstraints<BV, V, MI, C>;
  using Coefficient = typename Constraints::Coefficient;
  using ElementSeed = typename SecondaryBasis::GridView::Grid::template Codim<0>::EntitySeed;
  using GlobalCoordinate = typename SecondaryBasis::GridView::Grid::template Codim<0>::Geometry::GlobalCoordinate;
 
  static_assert(std::is_same_v<typename PrimaryBasis::MultiIndex, typename Constraints::MultiIndex>);
 
  // Interpolation weights smaller than the tolerance will be ignored
  const auto tol = 1e-5;
 
  auto localIsConstrained = std::vector<bool>();
  auto interpolationValues = std::vector<Coefficient>();
 
  auto primaryLocalView = primaryBasis.localView();
  auto secondaryLocalView = secondaryBasis.localView();
  auto secondaryFaceDOFs = Dune::Functions::subEntityDOFs(secondaryBasis);
 
  using Face = std::tuple<ElementSeed, int, GlobalCoordinate>;
 
  auto secondaryFaces = std::vector<Face>();
  for(const auto& secondaryIntersection : secondaryPatch)
    secondaryFaces.push_back(Face(
      secondaryIntersection.inside().seed(),
      secondaryIntersection.indexInInside(),
      secondaryIntersection.geometry().center()));
 
  for(const auto& primaryIntersection : primaryPatch)
  {
    auto primaryCenter = primaryIntersection.geometry().center();
    const auto& primaryElement = primaryIntersection.inside();
 
    for(const auto& secondaryFace : secondaryFaces)
    {
      const auto& [secondaryElementSeed, secondaryFaceIndex, secondaryCenter] = secondaryFace;
      const auto& secondaryElement = secondaryBasis.gridView().grid().entity(secondaryElementSeed);
 
      if ((primaryCenter - secondaryToPrimary(secondaryCenter)).two_norm() < tol)
      {
        primaryLocalView.bind(primaryElement);
        if (primaryLocalView.tree().size() == 0)
          break;
        secondaryLocalView.bind(secondaryElement);
        if (secondaryLocalView.tree().size() == 0)
          continue;
        secondaryFaceDOFs.bind(secondaryLocalView, secondaryFaceIndex, 1);
 
        const auto& primaryGeometry = primaryElement.geometry();
        const auto& secondaryGeometry = secondaryElement.geometry();
        auto localSecondaryToPrimary = [&](auto xLocalSecondary) {
          auto xGlobal = secondaryGeometry.global(xLocalSecondary);
          return primaryGeometry.local(secondaryToPrimary(xGlobal));
        };
 
        computeLocalInterpolationConstraints(primaryLocalView, secondaryLocalView, localSecondaryToPrimary, secondaryFaceDOFs, constraints);
 
        break;
      }
    }
  }
  constraints.resolveDependencies();
}
 
 
 
template<class BV, class V, class MI, class C, class Basis, class PrimaryPatch, class SecondaryPatch, class SecondaryToPrimary>
void computePeriodicConstraints(
  AffineConstraints<BV, V, MI, C>& constraints,
  const Basis& basis,
  const PrimaryPatch& primaryPatch,
  const SecondaryPatch& secondaryPatch,
  const SecondaryToPrimary& secondaryToPrimary)
{
  computePeriodicConstraints(constraints, basis, primaryPatch, basis, secondaryPatch, secondaryToPrimary);
}
 
 
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_CONSTRAINTS_PERIODICCONSTRAINTS_HH