GradientTest.cpp

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>

#include <amdis/AMDiS.hpp>

#include <amdis/AdaptiveGrid.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>

#include "Tests.hpp"

using namespace AMDiS;

template <std::size_t p>
void test()
{
  using HostGrid = Dune::YaspGrid<2>;
  using Basis1 = LagrangeBasis<HostGrid,p,p>;

  HostGrid hostGrid({1.0, 1.0}, {16,16});
  AdaptiveGrid<HostGrid> grid(hostGrid);
  auto const& gridView = grid.leafGridView();

  auto basis = Basis1::create(gridView);

  auto uVector = makeDOFVector(basis);

  auto u = valueOf(uVector);
  auto u_0 = valueOf(uVector, 0);
  auto u_1 = valueOf(uVector, 1);

  // eval a functor at global coordinates (at quadrature points)
  u_0 << evalAtQP([](auto const& x) { return x[0] + x[1]; });
  u_1 << evalAtQP([](auto const& x) { return -2*x[0] + 3*x[1]; });

  // test gradient
  AMDIS_TEST_APPROX((std::sqrt(integrate(unary_dot(gradientOf(u_0)), gridView))), std::sqrt(2.0));
  AMDIS_TEST_APPROX((std::sqrt(integrate(unary_dot(gradientOf(u_1)), gridView))), std::sqrt(13.0));

  // test divergence
  AMDIS_TEST_APPROX((integrate(divergenceOf(u), gridView)), 4.0);

  // test partial derivative
  double d0u_0 = integrate(partialDerivativeOf(u_0,0), gridView);
  double d1u_1 = integrate(partialDerivativeOf(u_1,1), gridView);
  AMDIS_TEST_APPROX(d0u_0, 1.0);
  AMDIS_TEST_APPROX(d1u_1, 3.0);

  auto gf_d0u_0 = makeGridFunction(partialDerivativeOf(u_0,0), gridView);
  auto lf_d0u_0 = localFunction(gf_d0u_0);
  for (auto const& e : elements(gridView)) {
    lf_d0u_0.bind(e);
    double v = lf_d0u_0(e.geometry().center());
    AMDIS_TEST_APPROX(v, 1.0);
    lf_d0u_0.unbind();
  }

  u_0 << evalAtQP([](auto const& x) { return std::sin(x[0]) + std::cos(x[1]); });
  u_1 << evalAtQP([](auto const& x) { return -std::sin(x[0]) * std::cos(x[1]); });

  AMDIS_TEST_APPROX((integrate(divergenceOf(u), gridView)),
                    (integrate(partialDerivativeOf(u_0,0) + partialDerivativeOf(u_1,1), gridView)));

  auto vVector(uVector);
  auto v = valueOf(vVector);
  auto v_0 = valueOf(vVector, 0);
  auto v_1 = valueOf(vVector, 1);

  // test gradient
  v << evalAtQP([](auto const& x) {
    return Dune::FieldVector<double,2>{std::cos(x[0]), -std::sin(x[1])};
  });
  AMDIS_TEST((std::sqrt(integrate(unary_dot(v - gradientOf(u_0)), gridView))) < 0.05);

  // test divergence
  v_0 << evalAtQP([](auto const& x) { return std::cos(x[0]) + std::sin(x[0])*std::sin(x[1]); });
  AMDIS_TEST((integrate(v_0 - divergenceOf(u), gridView)) < 0.05);

  // test partial derivative
  v_0 << evalAtQP([](auto const& x) { return std::cos(x[0]); });
  v_1 << evalAtQP([](auto const& x) { return std::sin(x[0])*std::sin(x[1]); });
  AMDIS_TEST((integrate(v_0 - partialDerivativeOf(u_0,0), gridView)) < 0.05);
  AMDIS_TEST((integrate(v_1 - partialDerivativeOf(u_1,1), gridView)) < 0.05);
}

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

  test<1u>();
  test<2u>();
  test<3u>();

  return report_errors();
}