p0p1interpolation.hh

Go to the documentation of this file.

// 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_FUNCTIONTOOLS_P0P1INTERPOLATION_HH
#define DUNE_FUFEM_FUNCTIONTOOLS_P0P1INTERPOLATION_HH
 
#include <vector>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/fufem/boundarypatch.hh>
 
template <class GridView, class Vector>
Vector interpolateP0ToP1(BoundaryPatch<GridView> const &patch,
                         Vector const &p0Data) {
  auto const &gridView = patch.gridView();
 
  std::vector<typename GridView::Grid::ctype> surroundingArea(
      gridView.size(GridView::dimension));
  std::fill(surroundingArea.begin(), surroundingArea.end(), 0);
 
  Vector ret(surroundingArea.size());
  ret = 0.0;
 
  auto const &indexSet = gridView.indexSet();
  for (auto const &intersection : patch) {
    auto const &inside = intersection.inside();
    auto refElement =
        Dune::ReferenceElements<double, GridView::dimension>::general(
            inside.type());
    auto const globalElementIndex = indexSet.index(inside);
 
    auto const localFaceIndex = intersection.indexInInside();
    auto const &face = inside.template subEntity<1>(localFaceIndex);
 
    auto const area = face.geometry().volume();
 
    int const numVertices =
        refElement.size(localFaceIndex, 1, GridView::dimension);
    for (int i = 0; i < numVertices; i++) {
      auto const localVertexIndex =
          refElement.subEntity(localFaceIndex, 1, i, GridView::dimension);
      auto const globalVertexIndex =
          indexSet.subIndex(inside, localVertexIndex, GridView::dimension);
 
      surroundingArea[globalVertexIndex] += area;
      ret[globalVertexIndex] += p0Data[globalElementIndex];
    }
  }
 
  for (size_t i = 0; i < ret.size(); ++i)
    if (patch.containsVertex(i))
      ret[i] /= surroundingArea[i];
 
  return ret;
}
 
template <class GridView, class Vector>
Vector interpolateP0ToP1(GridView const &gridView,
                         Vector const &p0Data) {
 
  unsigned int const dim = GridView::dimension;
  std::vector<typename GridView::Grid::ctype> surroundingArea(
      gridView.size(dim));
  std::fill(surroundingArea.begin(), surroundingArea.end(), 0);
 
  Vector ret(surroundingArea.size());
  ret = 0.0;
 
  auto const &indexSet = gridView.indexSet();
  for (auto const &element : elements(gridView)) {
    auto refElement =
        Dune::ReferenceElements<double, dim>::general(
            element.type());
    auto const globalElementIndex = indexSet.index(element);
    auto const area = element.geometry().volume();
 
    int const numVertices = refElement.size(dim);
    for (int i = 0; i < numVertices; i++) {
      auto const globalVertexIndex = indexSet.subIndex(element, i, dim);
 
      surroundingArea[globalVertexIndex] += area;
      ret[globalVertexIndex] += p0Data[globalElementIndex];
    }
  }
 
  for (size_t i = 0; i < ret.size(); ++i)
      ret[i] /= surroundingArea[i];
 
  return ret;
}
 
#endif