test1.cc 4.07 KB
Newer Older
1
2
3
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
4

5
#define TEST_EQ(expr1, expr2) test_exit( (expr1) == (expr2), "'" , #expr1 , "' != '" , #expr2 , "'" )
6
7
8
9
10
11
12
13
14
15

#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif

using namespace AMDiS;
namespace {
16

17
18
  using TestParam   = ProblemStatTraits<DIM, DOW, 1>;
  using TestProblem = ProblemStat<TestParam>;
19

20
21
22
23
  using MeshView = typename TestProblem::MeshView;
  using Geometry = typename MeshView::template Codim<0>::Geometry;
  using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;

24
25


26
27
  template <class FeSpace>
  class LocalEvaluation
28
  {
29
30
31
32
  public:
    LocalEvaluation(FeSpace const& feSpace)
      : localView(feSpace.localView())
    {}
33

34
35
36
37
38
    template <class Element>
    void bind(Element const& element)
    {
      localView.bind(element);
    }
39

40
41
42
43
44
45
46
    /// evaluate DOFVector at barycentric coordinates \p lambda in element that
    /// is bind to instance in \ref bind().
    template <class ValueType, class LocalCoord>
    ValueType eval(DOFVector<FeSpace, ValueType> const& vec, LocalCoord const& lambda)
    {
      auto const& localBasis = localView.tree().finiteElement().localBasis();
      localBasis.evaluateFunction(lambda, shapeValues);
47

48
      ValueType data = 0;
49
      for (std::size_t j = 0; j < shapeValues.size(); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
50
        const auto global_idx = localView.index(j);
51
52
        data += vec[global_idx] * shapeValues[j];
      }
53

54
55
      return data;
    }
56

57
58
  private:
    typename FeSpace::LocalView localView;
59

60
61
    std::vector<Dune::FieldVector<double,1> > shapeValues;
  };
62
63


64
  class AMDiSTester
65
  {
66
67
68
  public:
    // create a problemStat
    AMDiSTester(std::string name)
69
      : prob(make_shared<TestProblem>(name))
70
71
72
    {
      prob->initialize(INIT_ALL);
    }
73
74


75
76
77
    void test1() // CreateDOFVector
    {
      using FeSpace = typename TestProblem::template FeSpace<0>;
78

79
80
      // construction of DOFVectors
      DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
81
      auto vec2 = makeDOFVector(prob->getFeSpace(0_c), "vec2");
82

83
      // retrieving DOFVectors from problemStat
84
      auto vec3 = prob->getSolution(0_c);
85
86
      auto vec5 = prob->getSolution<0>();
      auto vec6 = prob->getSolution(index_<0>());
87

88
89
90
      auto vec7 = prob->getSolution()->getDOFVector(0_c);
      auto vec8 = prob->getSolution()->getDOFVector<0>();
      auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
91

92
93
94
      auto& sys_vec = *prob->getSolution();
      auto vec10 = sys_vec[0_c];
      auto vec11 = sys_vec[index_<0>()];
95

96
97
98
      // construction of SystemVector
      using FeSpaces = typename TestProblem::FeSpaces;
      SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
99
100
      auto sys_vec2 = makeSystemVector(*prob->getFeSpaces(), prob->getComponentNames());
      auto sys_vec3 = makeSystemVector(*prob->getFeSpaces(), "sys_vec");
101

102
103
104
      // retrieving systemVector from problemStat
      auto sys_vec4 = *prob->getSolution();
    }
105
106


107
108
109
110
    void test2() // FillDOFVector
    {
      auto& vec = prob->getSolution(0_c);
      auto const& feSpace = vec.getFeSpace();
111

112
113
      // interpolate function to dofvector
      vec.interpol([](auto const& x) { return x*x; });
114

115
116
117
      // test whether interpolation is fine
      using FeSpace = std::decay_t< decltype(feSpace) >;
      LocalEvaluation<FeSpace> evaluator(feSpace);
118

119
120
      for (auto const& element : elements(feSpace.gridView())) {
        evaluator.bind(element);
121

122
        auto const& refElement = RefElements::general(element.type());
123

124
125
        std::size_t nVertices = element.subEntities(DIM);
        for (std::size_t i = 0; i < nVertices; ++i) {
126
127
          auto pos = refElement.position(i, DIM);
          auto data = evaluator.eval(vec, pos);
128

129
130
131
          auto x = element.geometry().global(pos);
          TEST_EQ(data, x*x);
        }
132
133
      }
    }
134

135
  protected:
136
    shared_ptr<TestProblem> prob;
137
  };
138
139
140
141
142
143

} // end namespace

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

145
  auto tester = make_shared<AMDiSTester>("test");
146
147
  tester->test1();
  tester->test2();
148

149
  AMDiS::finalize();
150
  return 0;
151
}