DiscreteFunctionTest.cpp 9 KB
 1 ``````#include `````` 2 3 `````` #include `````` 4 5 ``````#include `````` 6 7 ``````#include `````` 8 9 ``````#include #include `````` 10 ``````#include `````` Praetorius, Simon committed Oct 23, 2018 11 ``````#include `````` Praetorius, Simon committed Feb 28, 2019 12 ``````#include `````` 13 14 15 16 17 `````` #include "Tests.hpp" using namespace AMDiS; `````` Müller, Felix committed Oct 09, 2019 18 ``````using ElliptParam = YaspGridBasis<2, 2, 2>; `````` 19 20 ``````using ElliptProblem = ProblemStat; `````` 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 ``````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 `````` 56 ``````template `````` 57 ``````bool compare(DOFVector const& U, DOFVector const& V) `````` 58 ``````{ `````` Müller, Felix committed Oct 09, 2019 59 60 61 62 `````` 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(); `````` 63 `````` if (&U.basis() != &V.basis()) `````` Müller, Felix committed Oct 09, 2019 64 `````` return false; `````` 65 66 `````` auto lv = basis.localView(); auto const& gv = basis.gridView(); `````` Müller, Felix committed Oct 09, 2019 67 68 69 70 71 72 `````` for (auto const& e : elements(gv)) { lv.bind(e); for (size_type i = 0; i < lv.size(); ++i) { auto multiIndex = lv.index(i); `````` Praetorius, Simon committed Dec 21, 2020 73 `````` if (!compare(U.get(multiIndex), V.get(multiIndex))) `````` Müller, Felix committed Oct 09, 2019 74 75 76 `````` return false; } } `````` 77 78 79 `````` return true; } `````` 80 81 82 ``````// compare discrete functions at quadrature points template bool compare(DiscreteFunction const& u, DiscreteFunction const& v) `````` 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 `````` 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)); } `````` 132 `````` `````` 133 134 135 136 `````` template void test2(Grid& grid, BasisFactory&& basis) { `````` 137 138 `````` using namespace Dune::Indices; `````` 139 `````` ProblemStat prob("test", grid, FWD(basis)); `````` 140 141 `````` prob.initialize(INIT_ALL); `````` Praetorius, Simon committed Mar 12, 2019 142 143 144 `````` auto U0 = *prob.solutionVector(); auto U1 = *prob.solutionVector(); auto U2 = *prob.solutionVector(); `````` 145 `````` `````` Praetorius, Simon committed Nov 05, 2020 146 147 148 `````` auto u0 = valueOf(U0); auto u1 = valueOf(U1); auto u2 = valueOf(U2); `````` Praetorius, Simon committed Dec 26, 2019 149 150 151 152 153 `````` // Test const and mutable construction auto const& U_c = *prob.solutionVector(); auto& U_m = *prob.solutionVector(); `````` 154 155 `````` DiscreteFunction u0_c{U_c.coefficients(), U_c.basis()}; DiscreteFunction u0_m{U_m.coefficients(), U_m.basis()}; `````` Praetorius, Simon committed Apr 14, 2020 156 `````` `````` Praetorius, Simon committed Nov 05, 2020 157 158 `````` auto u1_c = valueOf(U_c); auto u1_m = valueOf(U_m); `````` 159 `````` `````` Praetorius, Simon committed Dec 26, 2019 160 161 162 163 164 165 166 `````` // 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 167 168 169 170 171 172 173 174 175 `````` 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; `````` 176 `````` `````` 177 178 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 179 `````` `````` Müller, Felix committed Oct 09, 2019 180 181 182 183 `````` auto expr2 = [](auto const& x) { return Range{ 1 + 2*x[0] - 3*x[1], 2 + 2*x[0] - 3*x[1]}; }; `````` 184 185 186 187 188 `````` u0.interpolate_noalias(u2 + expr2); u1.interpolate(u1 + expr2); u2 += expr2; `````` 189 190 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 191 `````` `````` Müller, Felix committed Oct 09, 2019 192 193 194 195 `````` auto expr3 = [](auto const& x) { return Range{ 1 - 0.5*x[0] - 2*x[1], 2 - 0.5*x[0] - 2*x[1]}; }; `````` 196 197 198 199 200 `````` u0.interpolate_noalias(u2 - expr3); u1.interpolate(u1 - expr3); u2 -= expr3; `````` 201 202 `````` AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); `````` 203 `````` `````` Praetorius, Simon committed Nov 05, 2020 204 `````` auto du0 = derivativeOf(u0, tag::gradient{}); `````` 205 206 207 208 209 `````` auto lf0 = localFunction(u0); auto lf1 = localFunction(u1); auto dlf0 = derivativeOf(lf0, tag::gradient{}); auto dlf1 = derivative(lf1); `````` 210 211 `````` for (auto const& e : elements(prob.gridView())) { `````` 212 `````` lf0.bind(e); `````` 213 214 `````` auto geo = e.geometry(); auto local = referenceElement(geo).position(0,0); `````` 215 `````` auto y0 = lf0(local); `````` 216 217 `````` auto y1 = u0(geo.global(local)); `````` Müller, Felix committed Oct 09, 2019 218 `````` for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2) `````` 219 `````` AMDIS_TEST(compare(*it1, *it2)); `````` Müller, Felix committed Oct 09, 2019 220 `````` `````` 221 222 223 224 225 226 `````` // 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(); `````` 227 `````` `````` 228 229 230 231 232 233 `````` // should be bound independently to lf0 auto y1_copy = lf0_copy(local); lf0_copy.unbind(); dlf0.bind(e); auto g0 = dlf0(local); `````` 234 235 `````` auto g1 = du0(geo.global(local)); `````` Müller, Felix committed Oct 09, 2019 236 `````` for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2) `````` 237 `````` AMDIS_TEST(compare(*it1, *it2)); `````` 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 `````` // 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(); `````` 253 254 `````` } `````` Müller, Felix committed Oct 09, 2019 255 `````` auto V0 = makeDOFVector(prob.globalBasis()); `````` Praetorius, Simon committed Nov 05, 2020 256 `````` auto v0 = valueOf(V0); `````` Müller, Felix committed Oct 09, 2019 257 `````` v0 << expr1; `````` 258 `````` `````` Praetorius, Simon committed Nov 05, 2020 259 `````` // test DiscreteFunction `````` Praetorius, Simon committed Oct 14, 2020 260 261 `````` std::integral_constant _0; auto tp = makeTreePath(0); `````` Praetorius, Simon committed Mar 12, 2019 262 263 264 `````` auto W0 = *prob.solutionVector(); auto W1 = *prob.solutionVector(); auto W2 = *prob.solutionVector(); `````` Praetorius, Simon committed Nov 05, 2020 265 266 267 `````` auto w0 = valueOf(W0, 0); auto w1 = valueOf(W1, _0); auto w2 = valueOf(W2, tp); `````` 268 `````` `````` Praetorius, Simon committed Nov 05, 2020 269 `````` // test DiscreteFunction with (pre)treepath argument `````` Müller, Felix committed Oct 09, 2019 270 `````` auto expr = [](auto const& x) { return 1 + x[0] + x[1]; }; `````` Praetorius, Simon committed Mar 12, 2019 271 272 273 `````` auto W3 = *prob.solutionVector(); auto W4 = *prob.solutionVector(); auto W5 = *prob.solutionVector(); `````` Müller, Felix committed Oct 09, 2019 274 275 `````` auto W7 = *prob.solutionVector(); auto W8 = *prob.solutionVector(); `````` Praetorius, Simon committed Nov 05, 2020 276 277 278 `````` auto w3 = valueOf(W3, 0); auto w4 = valueOf(W4, _0); auto w5 = valueOf(W5, tp); `````` 279 `````` auto w6 = prob.solution(tp); `````` Praetorius, Simon committed Mar 12, 2019 280 `````` auto& W6 = *prob.solutionVector(); `````` 281 282 283 284 `````` w3 << expr; w4 << expr; w5 << expr; w6 << expr; `````` 285 286 287 `````` AMDIS_TEST( compare(W3, W4) ); AMDIS_TEST( compare(W3, W5) ); AMDIS_TEST( compare(W3, W6) ); `````` 288 `````` `````` Müller, Felix committed Oct 09, 2019 289 `````` // test interpolation on subbasis `````` Praetorius, Simon committed Nov 05, 2020 290 291 292 `````` auto w7 = valueOf(W7); auto w8_0 = valueOf(W8, 0); auto w8_1 = valueOf(W8, 1); `````` Müller, Felix committed Oct 09, 2019 293 294 295 `````` w7 << expr; w8_0 << expr; w8_1 << expr; `````` 296 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 `````` 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 346 `````` `````` 347 `````` return report_errors(); `````` 348 ``}``