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

5
#define TEST_EQ(expr1, expr2) AMDIS_TEST_EXIT( (expr1) == (expr2), "'" << #expr1 << "' != '" << #expr2 << "'" )
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
56
57
58
59
60
61
62
63
64
65
66

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

using namespace AMDiS;
namespace {
  
  using TestParam   = ProblemStatTraits<DIM, DOW, 1>;
  using TestProblem = ProblemStat<TestParam>;
  
  using MeshView = typename TestProblem::MeshView;
  using Geometry = typename MeshView::template Codim<0>::Geometry;
  using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;

  
  
  template <class FeSpace>
  class LocalEvaluation
  {   
  public:
    LocalEvaluation(FeSpace const& feSpace)
      : localView(feSpace.localView())
      , localIndexSet(feSpace.localIndexSet())
    {}
    
    template <class Element>
    void bind(Element const& element)
    {
      localView.bind(element);
      localIndexSet.bind(localView);
    }
    
    /// 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);
      
      ValueType data = 0;
      for (size_t j = 0; j < shapeValues.size(); ++j) {
        const auto global_idx = localIndexSet.index(j);
        data += vec[global_idx] * shapeValues[j];
      }
      
      return data;
    }
    
  private:
    typename FeSpace::LocalView localView;
    typename FeSpace::LocalIndexSet localIndexSet;
    
    std::vector<Dune::FieldVector<double,1> > shapeValues;
  };
  
  
67
  class AMDiSTester
68
  {
69
70
71
72
73
74
75
  public:
    // create a problemStat
    AMDiSTester(std::string name)
      : prob(std::make_shared<TestProblem>(name))
    {
      prob->initialize(INIT_ALL);
    }
76
77
    
  
78
79
80
81
82
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
    void test1() // CreateDOFVector
    {
      using FeSpace = typename TestProblem::template FeSpace<0>;
      
      // construction of DOFVectors
      DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
      auto vec2 = make_dofvector(prob->getFeSpace(0_c), "vec2");
      
      // retrieving DOFVectors from problemStat
      auto vec3 = prob->getSolution(0_c);    
      auto vec5 = prob->getSolution<0>();
      auto vec6 = prob->getSolution(index_<0>());
      
      auto vec7 = prob->getSolution()->getDOFVector(0_c);
      auto vec8 = prob->getSolution()->getDOFVector<0>();
      auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
      
      auto& sys_vec = *prob->getSolution();
      auto vec10 = sys_vec[0_c];
      auto vec11 = sys_vec[index_<0>()];
      
      // construction of SystemVector
      using FeSpaces = typename TestProblem::FeSpaces;
      SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
      auto sys_vec2 = make_systemvector(*prob->getFeSpaces(), prob->getComponentNames());
      auto sys_vec3 = make_systemvector(*prob->getFeSpaces(), "sys_vec");
      
      // retrieving systemVector from problemStat
      auto sys_vec4 = *prob->getSolution();
    }
108
  
109
110
111
112
113
  
    void test2() // FillDOFVector
    {
      auto& vec = prob->getSolution(0_c);
      auto const& feSpace = vec.getFeSpace();
114
      
115
116
      // interpolate function to dofvector
      vec.interpol([](auto const& x) { return x*x; });
117
      
118
119
120
121
122
123
124
125
      // test whether interpolation is fine
      using FeSpace = std::decay_t< decltype(feSpace) >;
      LocalEvaluation<FeSpace> evaluator(feSpace);
      
      for (auto const& element : elements(feSpace.gridView())) {
        evaluator.bind(element);
        
        auto const& refElement = RefElements::general(element.type());
126
        
127
128
129
130
131
132
133
134
        size_t nVertices = element.subEntities(DIM);
        for (size_t i = 0; i < nVertices; ++i) {
          auto pos = refElement.position(i, DIM);
          auto data = evaluator.eval(vec, pos);
          
          auto x = element.geometry().global(pos);
          TEST_EQ(data, x*x);
        }
135
136
      }
    }
137
138
139
140
    
  protected:
    std::shared_ptr<TestProblem> prob;
  };
141
142
143
144
145
146
147

} // end namespace

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);
  
148
149
150
  auto tester = std::make_shared<AMDiSTester>("test");
  tester->test1();
  tester->test2();
151
152
  
  AMDiS::finalize();
153
  return 0;
154
}