Dune Core Modules (2.8.0)

gridfunctiontest.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_TEST_GRIDFUNCTIONTEST_HH
4#define DUNE_FUNCTIONS_GRIDFUNCTIONS_TEST_GRIDFUNCTIONTEST_HH
5
7
8#include <dune/functions/common/functionconcepts.hh>
9#include <dune/functions/gridfunctions/gridfunction.hh>
10#include <dune/functions/gridfunctions/gridviewfunction.hh>
11#include <dune/functions/gridfunctions/gridviewentityset.hh>
12
13
14namespace Dune {
15namespace Functions {
16namespace Test {
17
18
19
20bool checkTrue(bool value, std::string error)
21{
22 if (not(value))
23 std::cout << "TEST FAILURE:" << error << std::endl;
24 return value;
25}
26
27
28
29template<class GridView, class F>
30double integrateGridViewFunction(const GridView& gridView, const F& f, unsigned int quadOrder)
31{
32 static const int dim = GridView::dimension;
33
34 double integral = 0;
35
36 auto fLocal = localFunction(f);
37
38 // Loop over elements and integrate over the function
39 for (const auto& e : elements(gridView))
40 {
41 auto geometry = e.geometry();
42
43 fLocal.bind(e);
44
45 // A quadrature rule
46 const auto& quad = QuadratureRules<double, dim>::rule(e.type(), quadOrder);
47
48 // Loop over all quadrature points
49 for ( size_t pt=0; pt < quad.size(); pt++ ) {
50
51 // Position of the current quadrature point in the reference element
52 auto quadPos = quad[pt].position();
53
54 // The multiplicative factor in the integral transformation formula
55 auto integrationElement = geometry.integrationElement(quadPos);
56
57 integral += fLocal(quadPos) * quad[pt].weight() * integrationElement;
58 }
59 fLocal.unbind();
60 }
61 return integral;
62}
63
64
65template<class GridView, class F>
66bool checkGridViewFunction(const GridView& gridView, const F& f, double exactIntegral, unsigned int quadOrder=1)
67{
68 bool passed = true;
69 double integral;
70
72
73 using EntitySet = typename F::EntitySet;
74 using Domain = typename EntitySet::GlobalCoordinate;
75 using Range = typename std::result_of<F(Domain)>::type;
76
77
78
79 std::cout << "Checking raw function f on grid view" << std::endl;
80 passed = checkTrue(isGridFunction<decltype(f), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
81 integral = integrateGridViewFunction(gridView, f, quadOrder);
82 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
83
84 std::cout << "Checking GridFunction<Range(Domain), EntitySet>(f) on grid view" << std::endl;
85 GridFunction<Range(Domain), EntitySet> f2 = f;
86 passed = checkTrue(isGridFunction<decltype(f2), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
87 integral = integrateGridViewFunction(gridView, f2, quadOrder);
88 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
89
90 std::cout << "Checking GridViewFunction<Range(Domain), GridView>(f) on grid view" << std::endl;
91 GridViewFunction<Range(Domain), GridView> f3 = f;
92 passed = checkTrue(isGridFunction<decltype(f3), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
93 integral = integrateGridViewFunction(gridView, f3, quadOrder);
94 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
95
96 std::cout << "Checking makeGridFunction(f) on grid view" << std::endl;
97 auto f4 = makeGridViewFunction(f, gridView);
98 passed = checkTrue(isGridFunction<decltype(f4), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
99 integral = integrateGridViewFunction(gridView, f4, quadOrder);
100 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
101
102 std::cout << "Checking GridViewFunction<Range(Domain), GridView>(makeGridFunction(f)) on grid view" << std::endl;
103 GridViewFunction<Range(Domain), GridView> f5 = f4;
104 integral = integrateGridViewFunction(gridView, f5, quadOrder);
105 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
106
107 std::cout << "Checking default constructed and assigned GridViewFunction<Range(Domain), GridView> on grid view" << std::endl;
108 GridViewFunction<Range(Domain), GridView> f6;
109 f6 = f5;
110 integral = integrateGridViewFunction(gridView, f6, quadOrder);
111 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
112
113 std::cout << "Checking reassigned GridViewFunction<Range(Domain), GridView> on grid view" << std::endl;
114 f6 = f3;
115 integral = integrateGridViewFunction(gridView, f6, quadOrder);
116 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
117
118 return passed;
119}
120
121
122
123
124} // namespace Test
125} // namespace Functions
126} // namespace Dune
127
128
129
130
131#endif // DUNE_FUNCTIONS_GRIDFUNCTIONS_TEST_GRIDFUNCTIONTEST_HH
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:280
static constexpr bool isGridFunction()
Check if F models the GridFunction concept with given signature and entity set.
Definition: functionconcepts.hh:268
void localFunction(DiscreteGlobalBasisFunction< TT... > &&t)=delete
Construction of local functions from a temporary DiscreteGlobalBasisFunction (forbidden)
@ dimension
The dimension of the grid.
Definition: gridview.hh:130
Dune namespace.
Definition: alignedallocator.hh:11
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)