DiscreteFunctionTest.cpp 4.91 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
#include <amdis/gridfunctions/DiscreteFunction.hpp>
10
#include <amdis/typetree/TreePath.hpp>
11
12
13
14
15

#include "Tests.hpp"

using namespace AMDiS;

16
using ElliptParam   = YaspGridBasis<2, 2, 2>;
17
18
19
20
21
using ElliptProblem = ProblemStat<ElliptParam>;

template <class GB, class T>
bool comp(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
{
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
  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 (basis.get() != V.basis().get())
    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 (std::abs(U.at(multiIndex) - V.at(multiIndex)) > AMDIS_TEST_TOL)
        return false;
    }
  }
40
41
42
43
44
  return true;
}

int main(int argc, char** argv)
{
45
  Environment env(argc, argv);
46
47
48
49
50
51
52

  using namespace Dune::Indices;

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

  // create 3 copies of the solution vector
53
54
55
  auto U0 = *prob.solutionVector();
  auto U1 = *prob.solutionVector();
  auto U2 = *prob.solutionVector();
56

57
58
59
60
61
62
63
64
65
66
67
68
69
70
  auto u0 = makeDiscreteFunction(U0);
  auto u1 = makeDiscreteFunction(U1);
  auto u2 = makeDiscreteFunction(U2);

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

#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  DiscreteFunction u0_c(U_c);
  DiscreteFunction u0_m(U_m);
#endif
  auto u1_c = makeDiscreteFunction(U_c);
  auto u1_m = makeDiscreteFunction(U_m);
71

72
73
74
75
76
77
78
  // 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);

79
80
81
82
83
84
85
86
87
88
  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;
89
90
91
92

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

93
94
95
96
  auto expr2 = [](auto const& x) {
    return Range{ 1 + 2*x[0] - 3*x[1],
                  2 + 2*x[0] - 3*x[1]};
  };
97
98
99
100
101
102
103
104

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

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

105
106
107
108
  auto expr3 = [](auto const& x) {
    return Range{ 1 - 0.5*x[0] - 2*x[1],
                  2 - 0.5*x[0] - 2*x[1]};
  };
109
110
111
112
113
114
115
116

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

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

117
  auto du0 = derivative(u0, tag::gradient{});
118
  auto localFct = localFunction(u0);
119
  auto dlocalFct = derivative(localFct, tag::gradient{});
120
121
122
123
124
125
126
127
  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));

128
129
130
    for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2)
      AMDIS_TEST(std::abs(*it1 - *it2) < 1.e-10);

131
132
133
134
135
136
    localFct.unbind();

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

137
138
    for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2)
      AMDIS_TEST(two_norm(*it1 - *it2) < 1.e-10);
139
140
141
    dlocalFct.unbind();
  }

142
  auto V0 = makeDOFVector<double>(prob.globalBasis());
143
  auto v0 = makeDiscreteFunction(V0);
144
  v0 << expr1;
145

146
  // test discreteFunction
147
148
149
  int preTP1 = 0;
  std::integral_constant<std::size_t, 0> preTP2;
  auto tp = treepath(preTP1);
150
151
152
  auto W0 = *prob.solutionVector();
  auto W1 = *prob.solutionVector();
  auto W2 = *prob.solutionVector();
153
154
155
156
  auto w0 = makeDiscreteFunction(W0, preTP1);
  auto w1 = makeDiscreteFunction(W1, preTP2);
  auto w2 = makeDiscreteFunction(W2, tp);

157
  // test discreteFunction with (pre)treepath argument
158
  auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
159
160
161
  auto W3 = *prob.solutionVector();
  auto W4 = *prob.solutionVector();
  auto W5 = *prob.solutionVector();
162
163
  auto W7 = *prob.solutionVector();
  auto W8 = *prob.solutionVector();
164
165
166
  auto w3 = makeDiscreteFunction(W3, preTP1);
  auto w4 = makeDiscreteFunction(W4, preTP2);
  auto w5 = makeDiscreteFunction(W5, tp);
167
  auto w6 = prob.solution(tp);
168
  auto& W6 = *prob.solutionVector();
169
170
171
172
173
174
175
176
  w3 << expr;
  w4 << expr;
  w5 << expr;
  w6 << expr;
  AMDIS_TEST( comp(W3, W4) );
  AMDIS_TEST( comp(W3, W5) );
  AMDIS_TEST( comp(W3, W6) );

177
  // test interpolation on subbasis
178
179
180
  auto w7 = makeDiscreteFunction(W7);
  auto w8_0 = makeDiscreteFunction(W8, 0);
  auto w8_1 = makeDiscreteFunction(W8, 1);
181
182
183
184
185
  w7 << expr;
  w8_0 << expr;
  w8_1 << expr;
  AMDIS_TEST( comp(W7, W8) );

186
187
  return 0;
}