1
2
3
4
5
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
67
68
69
70
71
72
73
74
75
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
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
|
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/functions/functionspacebases/taylorhoodbasis.hh>
#include <dune/functions/functionspacebases/test/basistest.hh>
using namespace Dune;
using namespace Dune::Functions;
int main (int argc, char* argv[]) try
{
Dune::MPIHelper::instance(argc, argv);
Dune::TestSuite test;
// Generate grid for testing
const int dim = 2;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
std::array<int,dim> elements = {{10, 10}};
GridType grid(l,elements);
typedef GridType::LeafGridView GridView;
typedef TaylorHoodBasis<GridView> Basis;
const GridView& gridView = grid.leafGridView();
Basis feBasis(gridView);
test.subTest(checkBasis(feBasis));
typedef Basis::MultiIndex MultiIndex;
// Sample the function f(x,y) = x on the grid vertices
// If we use that as the coefficients of a finite element function,
// we know its integral and can check whether quadrature returns
// the correct result
std::array<std::vector<double>,dim> x;
for (unsigned int i = 0; i<dim; i++)
x[i].resize(feBasis.size({i}));
// TODO: Implement interpolation properly using the global basis.
for (auto it = gridView.begin<dim>(); it != gridView.end<dim>(); ++it)
x[1][gridView.indexSet().index(*it)] = it->geometry().corner(0)[0];
// Objects required in the local context
auto localView = feBasis.localView();
auto localIndexSet = feBasis.localIndexSet();
std::vector<double> coefficients(localView.maxSize());
// Loop over elements and integrate over the function
double integral = 0;
for (auto it = gridView.begin<0>(); it != gridView.end<0>(); ++it)
{
localView.bind(*it);
localIndexSet.bind(localView);
// paranoia checks
assert(localView.size() == localIndexSet.size());
assert(&(localView.globalBasis()) == &(feBasis));
assert(&(localIndexSet.localView()) == &(localView));
// copy data from global vector
coefficients.resize(localIndexSet.size());
for (size_t i=0; i<localIndexSet.size(); i++)
{
MultiIndex mi = localIndexSet.index(i);
coefficients[i] = x[mi[0]][mi[1]];
}
// get access to the finite element
using namespace Dune::TypeTree::Indices;
typedef Basis::LocalView::Tree Tree;
auto& p_leaf = TypeTree::child(localView.tree(),_1);
auto& v_leaf = TypeTree::child(localView.tree(),_0,1);
auto& localFiniteElement = p_leaf.finiteElement();
// A quadrature rule
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(it->type(), 1);
// Loop over all quadrature points
for ( size_t pt=0; pt < quad.size(); pt++ ) {
// Position of the current quadrature point in the reference element
const FieldVector<double,dim>& quadPos = quad[pt].position();
// The multiplicative factor in the integral transformation formula
const double integrationElement = it->geometry().integrationElement(quadPos);
// Evaluate all shape function values at this point
std::vector<FieldVector<double,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
// Actually compute the vector entries
for (size_t i=0; i<localFiniteElement.localBasis().size(); i++)
{
integral += coefficients[p_leaf.localIndex(i)] * shapeFunctionValues[i] * quad[pt].weight() * integrationElement;
}
}
// unbind
localIndexSet.unbind();
localView.unbind();
}
std::cout << "Computed integral is " << integral << std::endl;
assert(std::abs(integral-0.5)< 1e-10);
return test.exit();
} catch ( Dune::Exception &e )
{
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch(...)
{
std::cerr << "Unknown exception thrown!" << std::endl;
return 1;
}
|
