DiscreteFunctionTest.cpp 9.04 KB
Newer Older
1
#include <config.h>
2
3

#include <iostream>
4
5
#include <type_traits>

6
7
#include <dune/common/version.hh>

8
9
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
10
11
#include <amdis/common/Apply.hpp>
#include <amdis/common/RecursiveForEach.hpp>
12
#include <amdis/gridfunctions/DiscreteFunction.hpp>
13
#include <amdis/typetree/TreePath.hpp>
14
15
16
17
18

#include "Tests.hpp"

using namespace AMDiS;

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

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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
57
template <class GB, class T>
58
bool compare(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
59
{
60
61
62
63
  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();
64
  if (&U.basis() != &V.basis())
65
    return false;
66
67
  auto lv = basis.localView();
  auto const& gv = basis.gridView();
68
69
70
71
72
73
  for (auto const& e : elements(gv))
  {
    lv.bind(e);
    for (size_type i = 0; i < lv.size(); ++i)
    {
      auto multiIndex = lv.index(i);
74
      if (!compare(U.at(multiIndex), V.at(multiIndex)))
75
76
77
        return false;
    }
  }
78
79
80
  return true;
}

81
82
83
// compare discrete functions at quadrature points
template <class... Args1, class... Args2>
bool compare(DiscreteFunction<Args1...> const& u, DiscreteFunction<Args2...> const& v)
84
{
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
  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));
}
133

134
135
136
137

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

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

143
144
145
  auto U0 = *prob.solutionVector();
  auto U1 = *prob.solutionVector();
  auto U2 = *prob.solutionVector();
146

147
148
149
  auto u0 = valueOf(U0);
  auto u1 = valueOf(U1);
  auto u2 = valueOf(U2);
150
151
152
153
154

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

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

158
159
  auto u1_c = valueOf(U_c);
  auto u1_m = valueOf(U_m);
160

161
162
163
164
165
166
167
  // 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);

168
169
170
171
172
173
174
175
176
  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;
177

178
179
  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );
180

181
182
183
184
  auto expr2 = [](auto const& x) {
    return Range{ 1 + 2*x[0] - 3*x[1],
                  2 + 2*x[0] - 3*x[1]};
  };
185
186
187
188
189

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

190
191
  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );
192

193
194
195
196
  auto expr3 = [](auto const& x) {
    return Range{ 1 - 0.5*x[0] - 2*x[1],
                  2 - 0.5*x[0] - 2*x[1]};
  };
197
198
199
200
201

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

202
203
  AMDIS_TEST( compare(U0, U1) );
  AMDIS_TEST( compare(U0, U2) );
204

205
  auto du0 = derivativeOf(u0, tag::gradient{});
206
207
208
209
210

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

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

222
223
224
225
226
227
    // 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();
228

229
230
231
232
233
234
    // should be bound independently to lf0
    auto y1_copy = lf0_copy(local);
    lf0_copy.unbind();

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

237
    for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2)
238
      AMDIS_TEST(compare(*it1, *it2));
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

    // 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();
254
255
  }

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

260
  // test DiscreteFunction
261
262
  std::integral_constant<std::size_t, 0> _0;
  auto tp = makeTreePath(0);
263
264
265
  auto W0 = *prob.solutionVector();
  auto W1 = *prob.solutionVector();
  auto W2 = *prob.solutionVector();
266
267
268
  auto w0 = valueOf(W0, 0);
  auto w1 = valueOf(W1, _0);
  auto w2 = valueOf(W2, tp);
269

270
  // test DiscreteFunction with (pre)treepath argument
271
  auto expr = [](auto const& x) { return 1 + x[0] + x[1]; };
272
273
274
  auto W3 = *prob.solutionVector();
  auto W4 = *prob.solutionVector();
  auto W5 = *prob.solutionVector();
275
276
  auto W7 = *prob.solutionVector();
  auto W8 = *prob.solutionVector();
277
278
279
  auto w3 = valueOf(W3, 0);
  auto w4 = valueOf(W4, _0);
  auto w5 = valueOf(W5, tp);
280
  auto w6 = prob.solution(tp);
281
  auto& W6 = *prob.solutionVector();
282
283
284
285
  w3 << expr;
  w4 << expr;
  w5 << expr;
  w6 << expr;
286
287
288
  AMDIS_TEST( compare(W3, W4) );
  AMDIS_TEST( compare(W3, W5) );
  AMDIS_TEST( compare(W3, W6) );
289

290
  // test interpolation on subbasis
291
292
293
  auto w7 = valueOf(W7);
  auto w8_0 = valueOf(W8, 0);
  auto w8_1 = valueOf(W8, 1);
294
295
296
  w7 << expr;
  w8_0 << expr;
  w8_1 << expr;
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
  AMDIS_TEST( compare(W7, W8) );
}

template <class GridView>
void test3(GridView const& gridView)
{
  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};
}


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());
347

348
  return report_errors();
349
}