DiscreteFunctionTest.cpp 9.04 KB
 1 ``````#include `````` 2 3 `````` #include `````` 4 5 ``````#include `````` 6 7 ``````#include `````` 8 9 ``````#include #include `````` 10 11 ``````#include #include `````` Praetorius, Simon committed Oct 23, 2018 12 ``````#include `````` Praetorius, Simon committed Feb 28, 2019 13 ``````#include `````` 14 15 16 17 18 `````` #include "Tests.hpp" using namespace AMDiS; `````` Müller, Felix committed Oct 09, 2019 19 ``````using ElliptParam = YaspGridBasis<2, 2, 2>; `````` 20 21 ``````using ElliptProblem = ProblemStat; `````` 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 ``````template , int> = 0> bool compare(T const& u, T const& v) { return std::abs(u - v) <= AMDIS_TEST_TOL * std::abs(u + v); } template bool compare(Dune::FieldVector const& u, Dune::FieldVector const& v) { for (int i = 0; i < n; ++i) if (!compare(u[i], v[i])) return false; return true; } template bool compare(Dune::FieldMatrix const& u, Dune::FieldMatrix const& v) { for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) if (!compare(u[i][j], v[i][j])) return false; return true; } template bool compare(Dune::TupleVector const& u, Dune::TupleVector const& v) { return Ranges::applyIndices([&](auto... i) { return (compare(u[i], v[i]) &&...); }); } // compare DOFVectors `````` 57 ``````template `````` 58 ``````bool compare(DOFVector const& U, DOFVector const& V) `````` 59 ``````{ `````` Müller, Felix committed Oct 09, 2019 60 61 62 63 `````` if (U.localSize() != V.localSize() || U.globalSize() != V.globalSize()) return false; using size_type = typename DOFVector::GlobalBasis::LocalView::size_type; auto const& basis = U.basis(); `````` 64 `````` if (&U.basis() != &V.basis()) `````` Müller, Felix committed Oct 09, 2019 65 `````` return false; `````` 66 67 `````` auto lv = basis.localView(); auto const& gv = basis.gridView(); `````` Müller, Felix committed Oct 09, 2019 68 69 70 71 72 73 `````` for (auto const& e : elements(gv)) { lv.bind(e); for (size_type i = 0; i < lv.size(); ++i) { auto multiIndex = lv.index(i); `````` 74 `````` if (!compare(U.at(multiIndex), V.at(multiIndex))) `````` Müller, Felix committed Oct 09, 2019 75 76 77 `````` return false; } } `````` 78 79 80 `````` return true; } `````` 81 82 83 ``````// compare discrete functions at quadrature points template bool compare(DiscreteFunction const& u, DiscreteFunction const& v) `````` 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 `````` auto uLocal = localFunction(u); auto vLocal = localFunction(v); auto const& gv = u.basis().gridView(); for (auto const& e : elements(gv)) { uLocal.bind(e); vLocal.bind(e); for (auto const& qp : Dune::QuadratureRules::rule(e.type(),4)) if (!compare(uLocal(qp.position()), vLocal(qp.position()))) return false; } return true; } template void test1(Grid& grid, BasisFactory&& basis) { ProblemStat prob("test", grid, FWD(basis)); prob.initialize(INIT_ALL); prob.solution() << [](auto const& x) { return x[0]*x[0] - x[1]*x[1] + 2*x[0]*x[1] + 4*x[0] - 3*x[1] + 2; }; // create copy of the solution vector auto U = *prob.solutionVector(); auto u = valueOf(U); auto v = prob.solution(); AMDIS_TEST(compare(u, v)); // Test const and mutable construction auto const& U_c = *prob.solutionVector(); auto& U_m = *prob.solutionVector(); DiscreteFunction u0_c{U_c.coefficients(), U_c.basis()}; DiscreteFunction u0_m{U_m.coefficients(), U_m.basis()}; auto u1_c = valueOf(U_c); auto u1_m = valueOf(U_m); AMDIS_TEST(compare(u0_c, u1_c)); AMDIS_TEST(compare(u0_m, u1_m)); // create a DiscreteFunction with a prescribed range type auto u2_c = valueOf(U_c); auto u3_c = u1_c.template child(); AMDIS_TEST(compare(u2_c, u3_c)); } `````` 133 `````` `````` 134 135 136 137 `````` template void test2(Grid& grid, BasisFactory&& basis) { `````` 138 139 `````` using namespace Dune::Indices; `````` 140 `````` ProblemStat prob("test", grid, FWD(basis)); `````` 141 142 `````` prob.initialize(INIT_ALL); `````` Praetorius, Simon committed Mar 12, 2019 143 144 145 `````` auto U0 = *prob.solutionVector(); auto U1 = *prob.solutionVector(); auto U2 = *prob.solutionVector(); `````` 146 `````` `````` Praetorius, Simon committed Nov 05, 2020 147 148 149 `````` auto u0 = valueOf(U0); auto u1 = valueOf(U1); auto u2 = valueOf(U2); `````` Praetorius, Simon committed Dec 26, 2019 150 151 152 153 154 `````` // Test const and mutable construction auto const& U_c = *prob.solutionVector(); auto& U_m = *prob.solutionVector(); `````` 155 156 `````` DiscreteFunction u0_c{U_c.coefficients(), U_c.basis()}; DiscreteFunction u0_m{U_m.coefficients(), U_m.basis()}; `````` Praetorius, Simon committed Apr 14, 2020 157 `````` `````` Praetorius, Simon committed Nov 05, 2020 158 159 `````` auto u1_c = valueOf(U_c); auto u1_m = valueOf(U_m); `````` 160 `````` `````` Praetorius, Simon committed Dec 26, 2019 161 162 163 164 165 166 167 `````` // sub-range view on DiscreteFunction auto su1_c = u1_c.child(); auto su1_c0 = u1_c.child(0); auto su1_m = u1_m.child(); auto su1_m0 = u1_m.child(0); `````` Müller, Felix committed Oct 09, 2019 168 169 170 171 172 173 174 175 176 `````` using Range = typename decltype(u0)::Range; auto expr1 = [](auto const& x) { return Range{ 1 + x[0] + x[1], 2 + x[0] + x[1]}; }; u0.interpolate_noalias(expr1); u1.interpolate(expr1); u2 << expr1; `````` 177 `````` `````` 178 179 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 180 `````` `````` Müller, Felix committed Oct 09, 2019 181 182 183 184 `````` auto expr2 = [](auto const& x) { return Range{ 1 + 2*x[0] - 3*x[1], 2 + 2*x[0] - 3*x[1]}; }; `````` 185 186 187 188 189 `````` u0.interpolate_noalias(u2 + expr2); u1.interpolate(u1 + expr2); u2 += expr2; `````` 190 191 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 192 `````` `````` Müller, Felix committed Oct 09, 2019 193 194 195 196 `````` auto expr3 = [](auto const& x) { return Range{ 1 - 0.5*x[0] - 2*x[1], 2 - 0.5*x[0] - 2*x[1]}; }; `````` 197 198 199 200 201 `````` u0.interpolate_noalias(u2 - expr3); u1.interpolate(u1 - expr3); u2 -= expr3; `````` 202 203 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 204 `````` `````` Praetorius, Simon committed Nov 05, 2020 205 `````` auto du0 = derivativeOf(u0, tag::gradient{}); `````` 206 207 208 209 210 `````` auto lf0 = localFunction(u0); auto lf1 = localFunction(u1); auto dlf0 = derivativeOf(lf0, tag::gradient{}); auto dlf1 = derivative(lf1); `````` 211 212 `````` for (auto const& e : elements(prob.gridView())) { `````` 213 `````` lf0.bind(e); `````` 214 215 `````` auto geo = e.geometry(); auto local = referenceElement(geo).position(0,0); `````` 216 `````` auto y0 = lf0(local); `````` 217 218 `````` auto y1 = u0(geo.global(local)); `````` Müller, Felix committed Oct 09, 2019 219 `````` for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2) `````` 220 `````` AMDIS_TEST(compare(*it1, *it2)); `````` Müller, Felix committed Oct 09, 2019 221 `````` `````` 222 223 224 225 226 227 `````` // create a copy of lf0 // NOTE: lf0_copy should be bound to the element already, since lf0 is bound auto lf0_copy = lf0; auto y0_copy = lf0_copy(local); lf0.unbind(); `````` 228 `````` `````` 229 230 231 232 233 234 `````` // should be bound independently to lf0 auto y1_copy = lf0_copy(local); lf0_copy.unbind(); dlf0.bind(e); auto g0 = dlf0(local); `````` 235 236 `````` auto g1 = du0(geo.global(local)); `````` Müller, Felix committed Oct 09, 2019 237 `````` for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2) `````` 238 `````` AMDIS_TEST(compare(*it1, *it2)); `````` 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 `````` // create a copy of dlf0 // NOTE: dlf0_copy should be bound to the element already, since dlf0 is bound auto dlf0_copy = dlf0; auto g0_copy = dlf0_copy(local); dlf0.unbind(); // should be bound independently to lf0 auto g1_copy = dlf0_copy(local); dlf0_copy.unbind(); dlf1.bind(e); auto g2 = dlf1(local); dlf1.unbind(); `````` 254 255 `````` } `````` Müller, Felix committed Oct 09, 2019 256 `````` auto V0 = makeDOFVector(prob.globalBasis()); `````` Praetorius, Simon committed Nov 05, 2020 257 `````` auto v0 = valueOf(V0); `````` Müller, Felix committed Oct 09, 2019 258 `````` v0 << expr1; `````` 259 `````` `````` Praetorius, Simon committed Nov 05, 2020 260 `````` // test DiscreteFunction `````` Praetorius, Simon committed Oct 14, 2020 261 262 `````` std::integral_constant _0; auto tp = makeTreePath(0); `````` Praetorius, Simon committed Mar 12, 2019 263 264 265 `````` auto W0 = *prob.solutionVector(); auto W1 = *prob.solutionVector(); auto W2 = *prob.solutionVector(); `````` Praetorius, Simon committed Nov 05, 2020 266 267 268 `````` auto w0 = valueOf(W0, 0); auto w1 = valueOf(W1, _0); auto w2 = valueOf(W2, tp); `````` 269 `````` `````` Praetorius, Simon committed Nov 05, 2020 270 `````` // test DiscreteFunction with (pre)treepath argument `````` Müller, Felix committed Oct 09, 2019 271 `````` auto expr = [](auto const& x) { return 1 + x[0] + x[1]; }; `````` Praetorius, Simon committed Mar 12, 2019 272 273 274 `````` auto W3 = *prob.solutionVector(); auto W4 = *prob.solutionVector(); auto W5 = *prob.solutionVector(); `````` Müller, Felix committed Oct 09, 2019 275 276 `````` auto W7 = *prob.solutionVector(); auto W8 = *prob.solutionVector(); `````` Praetorius, Simon committed Nov 05, 2020 277 278 279 `````` auto w3 = valueOf(W3, 0); auto w4 = valueOf(W4, _0); auto w5 = valueOf(W5, tp); `````` 280 `````` auto w6 = prob.solution(tp); `````` Praetorius, Simon committed Mar 12, 2019 281 `````` auto& W6 = *prob.solutionVector(); `````` 282 283 284 285 `````` w3 << expr; w4 << expr; w5 << expr; w6 << expr; `````` 286 287 288 `````` AMDIS_TEST( compare(W3, W4) ); AMDIS_TEST( compare(W3, W5) ); AMDIS_TEST( compare(W3, W6) ); `````` 289 `````` `````` Müller, Felix committed Oct 09, 2019 290 `````` // test interpolation on subbasis `````` Praetorius, Simon committed Nov 05, 2020 291 292 293 `````` auto w7 = valueOf(W7); auto w8_0 = valueOf(W8, 0); auto w8_1 = valueOf(W8, 1); `````` Müller, Felix committed Oct 09, 2019 294 295 296 `````` w7 << expr; w8_0 << expr; w8_1 << expr; `````` 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 `````` AMDIS_TEST( compare(W7, W8) ); } template void test3(GridView const& gridView) { using namespace Dune::Functions::BasisFactory; auto blockedBasis = makeBasis(gridView, power<2>(lagrange<2>())); std::vector> vec1(blockedBasis.size()); DiscreteFunction u1{vec1, blockedBasis}; Dune::BlockVector> vec2(blockedBasis.size()); DiscreteFunction u2{vec2, blockedBasis}; auto flatBasis = makeBasis(gridView, power<2>(lagrange<2>(), flatInterleaved())); std::vector vec3(flatBasis.size()); DiscreteFunction u3{vec3, flatBasis}; ISTLBlockVector vec4(flatBasis); DiscreteFunction u4{vec4, flatBasis}; } int main(int argc, char** argv) { Environment env(argc, argv); Dune::YaspGrid<2> grid({1.0, 1.0}, {2, 2}); using namespace Dune::Functions::BasisFactory; test1(grid, lagrange<1>()); test1(grid, composite(lagrange<1>(), lagrange<2>())); #if DUNE_VERSION_GT(DUNE_TYPETREE,2,6) test1(grid, lagrange(1)); test1(grid, power<2>(lagrange(2))); #endif test1>(grid, power<2>(lagrange<2>())); test1>(grid, power<2>(lagrange<2>())); #if DUNE_VERSION_GT(DUNE_TYPETREE,2,6) test2(grid, power<2>(lagrange(2))); #endif test2(grid, power<2>(lagrange<2>())); test3(grid.leafGridView()); `````` Müller, Felix committed Oct 09, 2019 347 `````` `````` 348 `````` return report_errors(); `````` 349 ``}``