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
|
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
/** @file
*
* \brief Test std::reference_wrapper<SomeContainer> as the CornerStorage in
* MultiLinearGeometry.
*/
#include <iostream>
#include <functional>
#include <limits>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
template<class ct>
struct TestGeometryTraits :
public Dune::MultiLinearGeometryTraits<ct>
{
template< int mydim, int cdim >
struct CornerStorage
{
typedef std::reference_wrapper<
const std::vector< Dune::FieldVector< ct, cdim > > > Type;
};
};
template<class ct, int mydim, int cdim>
using TestGeometry =
Dune::MultiLinearGeometry<ct, mydim, cdim, TestGeometryTraits<ct> >;
//! make a geometry that covers a kite shape
/**
* \return The returned geometry will reference a static array that hold the
* corner information. Because the order of destruction of static
* variables is a bit unclear, the returned value becomes invalid when
* main ends and shall not be used anymore.
*/
template<class ctype>
TestGeometry<ctype, 2, 2> kiteGeometry()
{
static const std::vector<Dune::FieldVector<ctype, 2> > kiteCorners{
{ 0, -1.5 }, // bottom
{ 1, 0 }, // right
{ -1, 0 }, // left
{ 0, 0.5 }, // top
};
return { Dune::ReferenceElements<ctype, 2>::cube(), kiteCorners };
}
template<class ctype, int dim>
bool expectCenter(const Dune::FieldVector<ctype, dim> &actual,
const Dune::FieldVector<ctype, dim> &expected)
{
bool match =
(actual - expected).two_norm2() < std::numeric_limits<ctype>::epsilon();
std::cout << (match ? "pass: " : "fail: ")
<< " Got: (" << actual << "), expected: (" << expected << ")"
<< std::endl;
return match;
}
int main()
{
bool pass = true;
std::cout << "making a geometry of a kite..." << std::endl;
auto geo = kiteGeometry<double>();
std::cout << "checking center of kite..." << std::endl;
pass &= expectCenter(geo.center(), { 0, -0.25 });
{
// the corners of an upward-pointing triangle
std::vector<Dune::FieldVector<double, 2> > triangleCorners{
{ -1, 0 }, // left
{ 1, 0 }, // right
{ 0, 1 }, // top
};
std::cout << "turn the geometry into an upward-pointing triangle..."
<< std::endl;
geo = {
Dune::ReferenceElements<double, 2>::simplex(),
triangleCorners,
};
std::cout << "checking center of upward-pointing triangle..." << std::endl;
pass &= expectCenter(geo.center(), { 0, 1.0/3 });
std::cout << "turning the geometry into a a leftward-pointing triangle "
<< "by moving the right corner to the bottom..." << std::endl
<< "(this is to show that the geometry really has a reference "
<< "to the coordinates, not a copy of them)" << std::endl;
triangleCorners[1] = { 0, -1 }; // move right corner to bottom
// NOTE: strictly speaking, the above modification invalidates geo.
// However, we're using MultiLinearGeometry here, which does cope with
// movement of vertices. We don't give any guarantee that it will
// continue to cope in the future, though. CachedMultiLinearGeometry
// already does not always handle this case, because it has no way of
// knowing that it would need to update its caches.
std::cout << "checking center of leftward-pointing triangle..."
<< std::endl;
pass &= expectCenter(geo.center(), { -1.0/3, 0 });
std::cout << "invalidating the geometry by letting the currently "
<< "referenced coordinate storage go out of scope..."
<< std::endl;
}
// triangleCorner is out of scope, so geo references freed storage and is
// now really invalid. The only remaining operations are assignment and
// destruction.
std::cout << "revalidating geometry by again assigning the geometry of a "
<< "kite..." << std::endl;
geo = kiteGeometry<double>();
std::cout << "checking center of kite (again)..." << std::endl;
pass &= expectCenter(geo.center(), { 0, -0.25 });
return pass ? 0 : 1;
}
|
