DUNE PDELab (git)

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
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_GRIDFUNCTIONS_TEST_GRIDFUNCTIONTEST_HH
8#define DUNE_FUNCTIONS_GRIDFUNCTIONS_TEST_GRIDFUNCTIONTEST_HH
9
10#include <type_traits>
11
13
14#include <dune/functions/common/functionconcepts.hh>
15#include <dune/functions/gridfunctions/gridfunction.hh>
16#include <dune/functions/gridfunctions/gridviewfunction.hh>
17#include <dune/functions/gridfunctions/gridviewentityset.hh>
18
19
20namespace Dune {
21namespace Functions {
22namespace Test {
23
24
25
26bool checkTrue(bool value, std::string error)
27{
28 if (not(value))
29 std::cout << "TEST FAILURE:" << error << std::endl;
30 return value;
31}
32
33
34
35template<class GridView, class F>
36double integrateGridViewFunction(const GridView& gridView, const F& f, unsigned int quadOrder)
37{
38 static const int dim = GridView::dimension;
39
40 double integral = 0;
41
42 auto fLocal = localFunction(f);
43
44 // Loop over elements and integrate over the function
45 for (const auto& e : elements(gridView))
46 {
47 auto geometry = e.geometry();
48
49 fLocal.bind(e);
50
51 if (not fLocal.bound())
52 DUNE_THROW(Exception, "Function is not bound after bind()");
53
54 // A quadrature rule
55 const auto& quad = QuadratureRules<double, dim>::rule(e.type(), quadOrder);
56
57 // Loop over all quadrature points
58 for ( size_t pt=0; pt < quad.size(); pt++ ) {
59
60 // Position of the current quadrature point in the reference element
61 auto quadPos = quad[pt].position();
62
63 // The multiplicative factor in the integral transformation formula
64 auto integrationElement = geometry.integrationElement(quadPos);
65
66 integral += fLocal(quadPos) * quad[pt].weight() * integrationElement;
67 }
68 fLocal.unbind();
69 }
70 return integral;
71}
72
73
74template<class GridView, class F>
75bool checkGridViewFunction(const GridView& gridView, const F& f, double exactIntegral, unsigned int quadOrder=1)
76{
77 bool passed = true;
78 double integral;
79
81
82 using EntitySet = typename F::EntitySet;
83 using Domain = typename EntitySet::GlobalCoordinate;
84 using Range = std::invoke_result_t<F, Domain>;
85
86
87
88 std::cout << "Checking raw function f on grid view" << std::endl;
89 passed = checkTrue(isGridFunction<decltype(f), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
90 integral = integrateGridViewFunction(gridView, f, quadOrder);
91 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
92
93 std::cout << "Checking GridFunction<Range(Domain), EntitySet>(f) on grid view" << std::endl;
94 GridFunction<Range(Domain), EntitySet> f2 = f;
95 passed = checkTrue(isGridFunction<decltype(f2), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
96 integral = integrateGridViewFunction(gridView, f2, quadOrder);
97 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
98
99 std::cout << "Checking GridViewFunction<Range(Domain), GridView>(f) on grid view" << std::endl;
100 GridViewFunction<Range(Domain), GridView> f3 = f;
101 passed = checkTrue(isGridFunction<decltype(f3), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
102 integral = integrateGridViewFunction(gridView, f3, quadOrder);
103 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
104
105 std::cout << "Checking makeGridFunction(f) on grid view" << std::endl;
106 auto f4 = makeGridViewFunction(f, gridView);
107 passed = checkTrue(isGridFunction<decltype(f4), Range(Domain), EntitySet>(), "Function does not model GridFunction concept");
108 integral = integrateGridViewFunction(gridView, f4, quadOrder);
109 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
110
111 std::cout << "Checking GridViewFunction<Range(Domain), GridView>(makeGridFunction(f)) on grid view" << std::endl;
112 GridViewFunction<Range(Domain), GridView> f5 = f4;
113 integral = integrateGridViewFunction(gridView, f5, quadOrder);
114 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
115
116 std::cout << "Checking default constructed and assigned GridViewFunction<Range(Domain), GridView> on grid view" << std::endl;
117 GridViewFunction<Range(Domain), GridView> f6;
118 f6 = f5;
119 integral = integrateGridViewFunction(gridView, f6, quadOrder);
120 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
121
122 std::cout << "Checking reassigned GridViewFunction<Range(Domain), GridView> on grid view" << std::endl;
123 f6 = f3;
124 integral = integrateGridViewFunction(gridView, f6, quadOrder);
125 passed = checkTrue(std::abs(integral-exactIntegral) < 1e-10, "Integral is "+std::to_string(integral)+" but should be "+std::to_string(exactIntegral));
126
127 return passed;
128}
129
130
131
132
133} // namespace Test
134} // namespace Functions
135} // namespace Dune
136
137
138
139
140#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:326
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
static constexpr bool isGridFunction()
Check if F models the GridFunction concept with given signature and entity set.
Definition: functionconcepts.hh:273
static constexpr int dimension
The dimension of the grid.
Definition: gridview.hh:134
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)