1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #include <config.h>
4
5 #include <vector>
6 #include <array>
7
8 #include <dune/common/parallel/mpihelper.hh>
9 #include <dune/common/reservedvector.hh>
10 #include <dune/common/classname.hh>
11 #include <dune/common/fvector.hh>
12 #include <dune/common/indices.hh>
13 #include <dune/common/tuplevector.hh>
14 #include <dune/common/test/testsuite.hh>
15
16 #include <dune/istl/bvector.hh>
17 #include <dune/istl/multitypeblockvector.hh>
18
19 #include <dune/functions/common/type_traits.hh>
20 #include <dune/functions/backends/concepts.hh>
21 #include <dune/functions/backends/istlvectorbackend.hh>
22
23
24 using namespace Dune;
25
26
27
28 /**
29 * \brief A Dummy global basis
30 *
31 * This is a mock class providing non-uniform size information.
32 * It's non-uniform in the sense that not all multi-indices
33 * have the same length.
34 */
35 template<std::size_t dim>
36 class GlobalBasisMoc
37 {
38 public:
39 using size_type = std::size_t;
40 using SizePrefix = Dune::ReservedVector<std::size_t, 3>;
41 using MultiIndex = Dune::ReservedVector<std::size_t, 3>;
42
43 /**
44 * \brief Construct from basis
45 */
GlobalBasisMoc()46 GlobalBasisMoc()
47 {}
48
49 /**
50 * \brief Return number possible values for next position in multi index
51 *
52 * This shall vanish. It's just here such that this can be used
53 * as size provider n place of the basis.
54 */
size(const SizePrefix & prefix) const55 size_type size(const SizePrefix& prefix) const
56 {
57 if (prefix.size() == 0)
58 return 2;
59 if (prefix.size() == 1)
60 {
61 if (prefix[0] == 0)
62 return 23;
63 if (prefix[0] == 1)
64 return 42;
65 }
66 if (prefix.size() == 2)
67 {
68 if (prefix[0] == 0)
69 return dim;
70 if (prefix[0] == 1)
71 return 0;
72 }
73 if (prefix.size() == 3)
74 return 0;
75 return 0;
76 }
77
operator size_type() const78 operator size_type () const
79 {
80 return 23*dim+42;
81 }
82
83 };
84
85
86 template<class Vector, class Coefficient, std::size_t dim, class MultiIndex>
checkISTLVectorBackend(std::string shortName="")87 Dune::TestSuite checkISTLVectorBackend(std::string shortName="")
88 {
89 Dune::TestSuite test(shortName);
90
91 using namespace Dune::Indices;
92 using Basis = GlobalBasisMoc<dim>;
93 using SizePrefix = typename Basis::SizePrefix;
94
95 Basis basis;
96
97 // Create raw vector
98 Vector x_raw;
99
100 // Create wrapped vector
101 auto x = Dune::Functions::istlVectorBackend(x_raw);
102
103 test.require(Dune::models<Dune::Functions::Concept::VectorBackend<Basis>, decltype(x)>(), "VectorBackend concept check")
104 << "Object returned by istlVectorBackend() does not model the VectorBackend concept";
105
106 // Resize wrapped vector using basis
107 x.resize(basis);
108
109 // Derive size information from vector
110 test.require(x_raw.size() == basis.size(SizePrefix{}), "resize check")
111 << "x_raw.size() is " << x_raw.size() << " but should be " << basis.size(SizePrefix{});
112
113 test.require(x_raw[_0].size() == basis.size(SizePrefix{0}), "resize check")
114 << "x_raw[_0].size() is " << x_raw[_0].size() << " but should be " << basis.size(SizePrefix{0});
115
116 for (std::size_t i=0; i<basis.size({0}); ++i)
117 test.require(x_raw[_0][i].size() == basis.size(SizePrefix{0,i}), "resize check")
118 << "x_raw[_0][" << i << "].size() is " << x_raw[_0][i].size() << " but should be " << basis.size(SizePrefix{0,i});
119
120 test.require(x_raw[_1].size() == basis.size(SizePrefix{1}), "resize check")
121 << "x_raw[_1].size() is " << x_raw[_0].size() << " but should be " << basis.size(SizePrefix{1});
122
123
124 // Assign values to each vector entry
125 for (std::size_t i=0; i<x_raw[_0].size(); ++i)
126 for (std::size_t j=0; j<x_raw[_0][i].size(); ++j)
127 x[MultiIndex{{0, i, j}}] = 0+i+j;
128 for (std::size_t i=0; i<x_raw[_1].size(); ++i)
129 x[MultiIndex{{1, i}}] = 1+i;
130
131
132 // Access vector entries via const reference
133 const auto& x_const = x;
134 for (std::size_t i=0; i<x_raw[_0].size(); ++i)
135 for (std::size_t j=0; j<x_raw[_0][i].size(); ++j)
136 {
137 test.check(x_const[MultiIndex{{0, i, j}}] == Coefficient(0+i+j))
138 << "x[{0," << i << "," << j << "}] contains wrong value";
139 }
140 for (std::size_t i=0; i<x_raw[_1].size(); ++i)
141 {
142 test.check(x_const[MultiIndex{{1, i}}] == Coefficient(1+i))
143 << "x[{1," << i << "}] contains wrong value";
144 }
145
146 return test;
147 }
148
149
main(int argc,char * argv[])150 int main (int argc, char *argv[]) try
151 {
152 // Set up MPI, if available
153 MPIHelper::instance(argc, argv);
154
155 ///////////////////////////////////
156 // Generate the grid
157 ///////////////////////////////////
158
159
160 Dune::TestSuite test;
161 {
162 using VelocityVector = std::vector<std::vector<double>>;
163 using PressureVector = std::vector<double>;
164 using Coefficient = double;
165 using Vector = Dune::TupleVector<VelocityVector, PressureVector>;
166 using MultiIndex = ReservedVector<std::size_t, 3>;
167 test.subTest(checkISTLVectorBackend<Vector, Coefficient, 2, MultiIndex>("TV<V<V<double>>, V<double>>"));
168 }
169
170 {
171 using VelocityVector = std::vector<Dune::BlockVector<Dune::FieldVector<double,1>>>;
172 using PressureVector = std::vector<Dune::FieldVector<double,1>>;
173 using Coefficient = double;
174 using Vector = Dune::TupleVector<VelocityVector, PressureVector>;
175 using MultiIndex = ReservedVector<std::size_t, 3>;
176 test.subTest(checkISTLVectorBackend<Vector, Coefficient, 2, MultiIndex>("TV<V<BV<FV<double,1>>>, V<FV<doule,1>>>"));
177 }
178
179 {
180 static const std::size_t dim = 5;
181 using VelocityVector = std::vector<std::array<Dune::FieldVector<double,1>,dim>>;
182 using PressureVector = std::vector<double>;
183 using Coefficient = double;
184 using Vector = Dune::TupleVector<VelocityVector, PressureVector>;
185 using MultiIndex = ReservedVector<std::size_t, 3>;
186 test.subTest(checkISTLVectorBackend<Vector, Coefficient, dim, MultiIndex>("TV<V<A<FV<double,1>,5>>, V<double>>"));
187 }
188
189 {
190 static const std::size_t dim = 5;
191 using VelocityVector = Dune::BlockVector<Dune::FieldVector<double,dim>>;
192 using PressureVector = Dune::BlockVector<Dune::FieldVector<double,1>>;
193 using Coefficient = double;
194 using Vector = Dune::MultiTypeBlockVector<VelocityVector, PressureVector>;
195 using MultiIndex = ReservedVector<std::size_t, 3>;
196 test.subTest(checkISTLVectorBackend<Vector, Coefficient, dim, MultiIndex>("MTBV<BV<FV<double,5>>, BV<FV<double,1>>>"));
197 }
198
199 {
200 static const std::size_t dim = 3;
201 using VelocityVector = std::vector<Dune::MultiTypeBlockVector<Dune::FieldVector<double,1>, double, Dune::FieldVector<double,1>>>;
202 using PressureVector = Dune::BlockVector<Dune::FieldVector<double,1>>;
203 using Coefficient = double;
204 using Vector = Dune::MultiTypeBlockVector<VelocityVector, PressureVector>;
205 using MultiIndex = ReservedVector<std::size_t, 3>;
206 test.subTest(checkISTLVectorBackend<Vector, Coefficient, dim, MultiIndex>("MTBV<V<MTBV<FV<double,1>, double, FV<double,1>>>, BV<FV<double,1>>"));
207 }
208
209
210 return test.exit();
211 }
212 // Error handling
213 catch (Exception& e) {
214 std::cout << e.what() << std::endl;
215 return 1;
216 }
217