DiscreteFunctionTest.cpp

#include <config.h>

#include <iostream>
#include <type_traits>

#include <dune/common/version.hh>

#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Apply.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/typetree/TreePath.hpp>

#include "Tests.hpp"

using namespace AMDiS;

using ElliptParam   = YaspGridBasis<2, 2, 2>;
using ElliptProblem = ProblemStat<ElliptParam>;

template <class T,
  std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
bool compare(T const& u, T const& v)
{
  return std::abs(u - v) <= AMDIS_TEST_TOL * std::abs(u + v);
}

template <class T, int n>
bool compare(Dune::FieldVector<T,n> const& u, Dune::FieldVector<T,n> const& v)
{
  for (int i = 0; i < n; ++i)
    if (!compare(u[i], v[i]))
      return false;
  return true;
}

template <class T, int n, int m>
bool compare(Dune::FieldMatrix<T,n,m> const& u, Dune::FieldMatrix<T,n,m> const& v)
{
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < m; ++j)
      if (!compare(u[i][j], v[i][j]))
        return false;
  return true;
}

template <class... T>
bool compare(Dune::TupleVector<T...> const& u, Dune::TupleVector<T...> const& v)
{
  return Ranges::applyIndices<sizeof...(T)>([&](auto... i) {
    return (compare(u[i], v[i]) &&...);
  });
}

// compare DOFVectors
template <class GB, class T>
bool compare(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
{
  if (U.localSize() != V.localSize() || U.globalSize() != V.globalSize())
    return false;
  using size_type = typename DOFVector<GB,T>::GlobalBasis::LocalView::size_type;
  auto const& basis = U.basis();
  if (&U.basis() != &V.basis())
    return false;
  auto lv = basis.localView();
  auto const& gv = basis.gridView();
  for (auto const& e : elements(gv))
  {
    lv.bind(e);
    for (size_type i = 0; i < lv.size(); ++i)
    {
      auto multiIndex = lv.index(i);
      if (!compare(U.get(multiIndex), V.get(multiIndex)))
        return false;
    }
  }
  return true;
}

// compare discrete functions at quadrature points
template <class... Args1, class... Args2>
bool compare(DiscreteFunction<Args1...> const& u, DiscreteFunction<Args2...> const& v)
{
  auto uLocal = localFunction(u);
  auto vLocal = localFunction(v);

  auto const& gv = u.basis().gridView();
  for (auto const& e : elements(gv)) {
    uLocal.bind(e);
    vLocal.bind(e);
    for (auto const& qp : Dune::QuadratureRules<double,TYPEOF(gv)::dimension>::rule(e.type(),4))
      if (!compare(uLocal(qp.position()), vLocal(qp.position())))
        return false;
  }
  return true;
}


template <class Range = void, class Grid, class BasisFactory>
void test1(Grid& grid, BasisFactory&& basis)
{
  ProblemStat prob("test", grid, FWD(basis));
  prob.initialize(INIT_ALL);

  prob.solution() << [](auto const& x) { return x[0]*x[0] - x[1]*x[1] + 2*x[0]*x[1] + 4*x[0] - 3*x[1] + 2; };

  // create copy of the solution vector
  auto U = *prob.solutionVector();
  auto u = valueOf(U);
  auto v = prob.solution();

  AMDIS_TEST(compare(u, v));

  // Test const and mutable construction
  auto const& U_c = *prob.solutionVector();
  auto& U_m = *prob.solutionVector();

  DiscreteFunction u0_c{U_c.coefficients(), U_c.basis()};
  DiscreteFunction u0_m{U_m.coefficients(), U_m.basis()};

  auto u1_c = valueOf(U_c);
  auto u1_m = valueOf(U_m);

  AMDIS_TEST(compare(u0_c, u1_c));
  AMDIS_TEST(compare(u0_m, u1_m));

  // create a DiscreteFunction with a prescribed range type
  auto u2_c = valueOf<Range>(U_c);
  auto u3_c = u1_c.template child<Range>();
  AMDIS_TEST(compare(u2_c, u3_c));
}


template <class Grid, class BasisFactory>
void test2(Grid& grid, BasisFactory&& basis)
{
  using namespace Dune::Indices;

  ProblemStat prob("test", grid, FWD(basis));
  prob.initialize(INIT_ALL);

  auto U0 = *prob.solutionVector();
  auto U1 = *prob.solutionVector();
  auto U2 = *prob.solutionVector();

  auto u0 = valueOf(U0);
  auto u1 = valueOf(U1);
  auto u2 = valueOf(U2);

  // Test const and mutable construction
  auto const& U_c = *prob.solutionVector();
  auto& U_m = *prob.solutionVector();

  DiscreteFunction u0_c{U_c.coefficients(), U_c.basis()};
  DiscreteFunction u0_m{U_m.coefficients(), U_m.basis()};

  auto u1_c = valueOf(U_c);
  auto u1_m = valueOf(U_m);

  // sub-range view on DiscreteFunction
  auto su1_c = u1_c.child();
  auto su1_c0 = u1_c.child(0);

  auto su1_m = u1_m.child();
  auto su1_m0 = u1_m.child(0);

  using Range = typename decltype(u0)::Range;
  auto expr1 = [](auto const& x) {
    return Range{ 1 + x[0] + x[1],
                  2 + x[0] + x[1]};
  };

  u0.interpolate_noalias(expr1);
  u1.interpolate(expr1);
  u2 << expr1;

  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );

  auto expr2 = [](auto const& x) {
    return Range{ 1 + 2*x[0] - 3*x[1],
                  2 + 2*x[0] - 3*x[1]};
  };

  u0.interpolate_noalias(u2 + expr2);
  u1.interpolate(u1 + expr2);
  u2 += expr2;

  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );

  auto expr3 = [](auto const& x) {
    return Range{ 1 - 0.5*x[0] - 2*x[1],
                  2 - 0.5*x[0] - 2*x[1]};
  };

  u0.interpolate_noalias(u2 - expr3);
  u1.interpolate(u1 - expr3);
  u2 -= expr3;

  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );

  auto du0 = derivativeOf(u0, tag::gradient{});

  auto lf0 = localFunction(u0);
  auto lf1 = localFunction(u1);
  auto dlf0 = derivativeOf(lf0, tag::gradient{});
  auto dlf1 = derivative(lf1);
  for (auto const& e : elements(prob.gridView()))
  {
    lf0.bind(e);
    auto geo = e.geometry();
    auto local = referenceElement(geo).position(0,0);
    auto y0 = lf0(local);
    auto y1 = u0(geo.global(local));

    for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2)
      AMDIS_TEST(compare(*it1, *it2));

    // create a copy of lf0
    // NOTE: lf0_copy should be bound to the element already, since lf0 is bound
    auto lf0_copy = lf0;
    auto y0_copy = lf0_copy(local);

    lf0.unbind();

    // should be bound independently to lf0
    auto y1_copy = lf0_copy(local);
    lf0_copy.unbind();

    dlf0.bind(e);
    auto g0 = dlf0(local);
    auto g1 = du0(geo.global(local));

    for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2)
      AMDIS_TEST(compare(*it1, *it2));

    // create a copy of dlf0
    // NOTE: dlf0_copy should be bound to the element already, since dlf0 is bound
    auto dlf0_copy = dlf0;
    auto g0_copy = dlf0_copy(local);

    dlf0.unbind();

    // should be bound independently to lf0
    auto g1_copy = dlf0_copy(local);
    dlf0_copy.unbind();

    dlf1.bind(e);
    auto g2 = dlf1(local);
    dlf1.unbind();
  }

  auto V0 = makeDOFVector<double>(prob.globalBasis());
  auto v0 = valueOf(V0);
  v0 << expr1;

  // test DiscreteFunction
  std::integral_constant<std::size_t, 0> _0;
  auto tp = makeTreePath(0);
  auto W0 = *prob.solutionVector();
  auto W1 = *prob.solutionVector();
  auto W2 = *prob.solutionVector();
  auto w0 = valueOf(W0, 0);
  auto w1 = valueOf(W1, _0);
  auto w2 = valueOf(W2, tp);

  // test DiscreteFunction with (pre)treepath argument
  auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
  auto W3 = *prob.solutionVector();
  auto W4 = *prob.solutionVector();
  auto W5 = *prob.solutionVector();
  auto W7 = *prob.solutionVector();
  auto W8 = *prob.solutionVector();
  auto w3 = valueOf(W3, 0);
  auto w4 = valueOf(W4, _0);
  auto w5 = valueOf(W5, tp);
  auto w6 = prob.solution(tp);
  auto& W6 = *prob.solutionVector();
  w3 << expr;
  w4 << expr;
  w5 << expr;
  w6 << expr;
  AMDIS_TEST( compare(W3, W4) );
  AMDIS_TEST( compare(W3, W5) );
  AMDIS_TEST( compare(W3, W6) );

  // test interpolation on subbasis
  auto w7 = valueOf(W7);
  auto w8_0 = valueOf(W8, 0);
  auto w8_1 = valueOf(W8, 1);
  w7 << expr;
  w8_0 << expr;
  w8_1 << expr;
  AMDIS_TEST( compare(W7, W8) );
}

template <class GridView>
void test3(GridView const& gridView)
{
#if AMDIS_BACKEND == AMDIS_BACKEND_ISTL
  using namespace Dune::Functions::BasisFactory;
  auto blockedBasis = makeBasis(gridView, power<2>(lagrange<2>()));

  std::vector<Dune::FieldVector<double,2>> vec1(blockedBasis.size());
  DiscreteFunction u1{vec1, blockedBasis};

  Dune::BlockVector<Dune::FieldVector<double,2>> vec2(blockedBasis.size());
  DiscreteFunction u2{vec2, blockedBasis};

  auto flatBasis = makeBasis(gridView, power<2>(lagrange<2>(), flatInterleaved()));

  std::vector<double> vec3(flatBasis.size());
  DiscreteFunction u3{vec3, flatBasis};

  ISTLBlockVector<double> vec4(flatBasis);
  DiscreteFunction u4{vec4, flatBasis};
#endif
}


int main(int argc, char** argv)
{
  Environment env(argc, argv);

  Dune::YaspGrid<2> grid({1.0, 1.0}, {2, 2});

  using namespace Dune::Functions::BasisFactory;
  test1(grid, lagrange<1>());
  test1(grid, composite(lagrange<1>(), lagrange<2>()));

#if DUNE_VERSION_GT(DUNE_TYPETREE,2,6)
  test1(grid, lagrange(1));
  test1(grid, power<2>(lagrange(2)));
#endif

  test1<Dune::FieldVector<double,2>>(grid, power<2>(lagrange<2>()));
  test1<Dune::FieldVector<float,2>>(grid, power<2>(lagrange<2>()));

#if DUNE_VERSION_GT(DUNE_TYPETREE,2,6)
  test2(grid, power<2>(lagrange(2)));
#endif

  test2(grid, power<2>(lagrange<2>()));

  test3(grid.leafGridView());

  return report_errors();
}