DiscreteFunctionTest.cpp 9.06 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
#include <amdis/common/Apply.hpp>
11
#include <amdis/gridfunctions/DiscreteFunction.hpp>
12
#include <amdis/typetree/TreePath.hpp>
13
14
15
16
17

#include "Tests.hpp"

using namespace AMDiS;

18
using ElliptParam   = LagrangeBasis<Dune::YaspGrid<2>,1,1>;
19
20
using ElliptProblem = ProblemStat<ElliptParam>;

21
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
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
56
template <class GB, class T>
57
bool compare(DOFVector<GB,T> const& U, DOFVector<GB,T> const& V)
58
{
59
60
61
62
  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();
63
  if (&U.basis() != &V.basis())
64
    return false;
65
66
  auto lv = basis.localView();
  auto const& gv = basis.gridView();
67
68
69
70
71
72
  for (auto const& e : elements(gv))
  {
    lv.bind(e);
    for (size_type i = 0; i < lv.size(); ++i)
    {
      auto multiIndex = lv.index(i);
73
      if (!compare(U.get(multiIndex), V.get(multiIndex)))
74
75
76
        return false;
    }
  }
77
78
79
  return true;
}

80
81
82
// compare discrete functions at quadrature points
template <class... Args1, class... Args2>
bool compare(DiscreteFunction<Args1...> const& u, DiscreteFunction<Args2...> const& v)
83
{
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
  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));
}
132

133
134
135
136

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

289
  // test interpolation on subbasis
290
291
292
  auto w7 = valueOf(W7);
  auto w8_0 = valueOf(W8, 0);
  auto w8_1 = valueOf(W8, 1);
293
294
295
  w7 << expr;
  w8_0 << expr;
  w8_1 << expr;
296
297
298
299
300
301
  AMDIS_TEST( compare(W7, W8) );
}

template <class GridView>
void test3(GridView const& gridView)
{
302
#if AMDIS_BACKEND == AMDIS_BACKEND_ISTL
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
  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};
319
#endif
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
347
}


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

349
  return report_errors();
350
}