DiscreteFunctionTest.cpp 3.32 KB
Newer Older
1 2 3 4
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>
5 6
#include <type_traits>

7 8
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
9 10
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>
11
#include <amdis/typetree/TreePath.hpp>
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

#include "Tests.hpp"

using namespace AMDiS;

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

template <class GB, class T>
bool comp(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
{
  if (U.size() != V.size())
    return false;

  using Int = typename DOFVector<GB,T>::size_type;
  for (Int i = 0; i < U.size(); ++i) {
    if (U[i] != V[i])
      return false;
  }
  return true;
}

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

  using namespace Dune::Indices;

  ElliptProblem prob("ellipt");
  prob.initialize(INIT_ALL);

  // create 3 copies of the solution vector
44 45 46
  auto U0 = *prob.solutionVector();
  auto U1 = *prob.solutionVector();
  auto U2 = *prob.solutionVector();
47 48 49 50 51

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

52
  auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
53 54 55 56 57 58 59
  u0.interpolate_noalias(expr);
  u1.interpolate(expr);
  u2 << expr;

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

60
  auto expr2 = [](auto const& x) { return 1 + 2*x[0] - 3*x[1]; };
61 62 63 64 65 66 67 68

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

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

69
  auto expr3 = [](auto const& x) { return 1 - 0.5*x[0] - 2*x[1]; };
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

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

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

  auto du0 = derivative(u0);
  auto localFct = localFunction(u0);
  auto dlocalFct = derivative(localFct);
  for (auto const& e : elements(prob.gridView()))
  {
    localFct.bind(e);
    auto geo = e.geometry();
    auto local = referenceElement(geo).position(0,0);
    auto y0 = localFct(local);
    auto y1 = u0(geo.global(local));

    AMDIS_TEST(std::abs(y0 - y1) < 1.e-10);
    localFct.unbind();

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

    AMDIS_TEST(two_norm(g0 - g1) < 1.e-10);
    dlocalFct.unbind();
  }

100
  auto V0 = makeDOFVector<double>(*prob.globalBasis());
101 102 103
  auto v0 = makeDOFVectorView(V0);
  v0 << expr;

104 105 106 107
  // test makeDiscreteFunction
  int preTP1 = 0;
  std::integral_constant<std::size_t, 0> preTP2;
  auto tp = treepath(preTP1);
108 109 110
  auto W0 = *prob.solutionVector();
  auto W1 = *prob.solutionVector();
  auto W2 = *prob.solutionVector();
111 112 113 114 115
  auto w0 = makeDiscreteFunction(W0, preTP1);
  auto w1 = makeDiscreteFunction(W1, preTP2);
  auto w2 = makeDiscreteFunction(W2, tp);

  // test makeDOFVectorView with (pre)treepath argument
116 117 118
  auto W3 = *prob.solutionVector();
  auto W4 = *prob.solutionVector();
  auto W5 = *prob.solutionVector();
119 120 121 122
  auto w3 = makeDOFVectorView(W3, preTP1);
  auto w4 = makeDOFVectorView(W4, preTP2);
  auto w5 = makeDOFVectorView(W5, tp);
  auto w6 = prob.solution(tp);
123
  auto& W6 = *prob.solutionVector();
124 125 126 127 128 129 130 131
  w3 << expr;
  w4 << expr;
  w5 << expr;
  w6 << expr;
  AMDIS_TEST( comp(W3, W4) );
  AMDIS_TEST( comp(W3, W5) );
  AMDIS_TEST( comp(W3, W6) );

132 133 134
  AMDiS::finalize();
  return 0;
}