GradientTest.cpp 4.24 KB
Newer Older
1
2
3
4
5
6
// -*- 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>
7
8

#include <amdis/AdaptiveGrid.hpp>
9
10
11
#include <amdis/LinearAlgebra.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/Integrate.hpp>
12
#include <amdis/gridfunctions/DiscreteFunction.hpp>
13

14
15
16
17
18
19
20
21
#if HAVE_DUNE_UGGRID
  #include <dune/grid/uggrid.hh>
#else
  #include <dune/grid/yaspgrid.hh>
#endif

#include <dune/grid/utility/structuredgridfactory.hh>

22
23
24
25
#include "Tests.hpp"

using namespace AMDiS;

26
27
template <std::size_t p, class HostGrid>
void test(HostGrid& hostGrid)
28
{
29
  std::cout << "[p = " << p << "]" << std::endl;
30
  using Basis1 = LagrangeBasis<HostGrid,p,p>;
Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
  AdaptiveGrid<HostGrid> grid(hostGrid);
  auto const& gridView = grid.leafGridView();
33
34
35
36
37

  auto basis = Basis1::create(gridView);

  auto uVector = makeDOFVector(basis);

38
39
40
  auto u = valueOf(uVector);
  auto u_0 = valueOf(uVector, 0);
  auto u_1 = valueOf(uVector, 1);
41
42
43
44
45
46

  // 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
47
48
  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));
49
50

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

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

59
  auto gf_d0u_0 = makeGridFunction(partialDerivativeOf(u_0,0), gridView);
60
61
62
63
64
65
66
67
68
69
70
  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]); });

71
72
  AMDIS_TEST_APPROX((integrate(divergenceOf(u), gridView)),
                    (integrate(partialDerivativeOf(u_0,0) + partialDerivativeOf(u_1,1), gridView)));
73
74

  auto vVector(uVector);
75
76
77
  auto v = valueOf(vVector);
  auto v_0 = valueOf(vVector, 0);
  auto v_1 = valueOf(vVector, 1);
78
79

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

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

  // 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]); });
92
93
  AMDIS_TEST((integrate(v_0 - partialDerivativeOf(u_0,0), gridView)) < 0.05);
  AMDIS_TEST((integrate(v_1 - partialDerivativeOf(u_1,1), gridView)) < 0.05);
94
95
96
97
98
99
}

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

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
133
134
135
#if HAVE_DUNE_UGGRID
  std::cout << std::endl << "UGGrid<Simplex>..." << std::endl;
  {
    using HostGrid = Dune::UGGrid<2>;
    auto numElements = Dune::filledArray<2,unsigned int>(16);
    auto hostGridPtr = Dune::StructuredGridFactory<HostGrid>::createSimplexGrid({0.0,0.0}, {1.0,1.0}, numElements);

    test<1u>(*hostGridPtr);
    test<2u>(*hostGridPtr);
    test<3u>(*hostGridPtr);
    test<4u>(*hostGridPtr);
  }

  std::cout << std::endl << "UGGrid<Cube>..." << std::endl;
  {
    using HostGrid = Dune::UGGrid<2>;
    auto numElements = Dune::filledArray<2,unsigned int>(16);
    auto hostGridPtr = Dune::StructuredGridFactory<HostGrid>::createCubeGrid({0.0,0.0}, {1.0,1.0}, numElements);

    test<1u>(*hostGridPtr);
    test<2u>(*hostGridPtr);
    test<3u>(*hostGridPtr);
    test<4u>(*hostGridPtr);
  }
#else
  std::cout << "YaspGrid..." << std::endl;
  {
    using HostGrid = Dune::YaspGrid<2>;
    HostGrid hostGrid({1.0, 1.0}, {16,16});

    test<1u>(hostGrid);
    test<2u>(hostGrid);
    test<3u>(hostGrid);
    test<4u>(hostGrid);
  }
#endif
136
137
138

  return report_errors();
}