#include #include #include #include #include #include #include #include #include #include #include "Tests.hpp" using namespace AMDiS; using ElliptParam = YaspGridBasis<2, 2, 2>; using ElliptProblem = ProblemStat; 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 template bool compare(DOFVector const& U, DOFVector const& V) { 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(); if (&U.basis() != &V.basis()) return false; auto lv = basis.localView(); auto const& gv = basis.gridView(); for (auto const& e : elements(gv)) { lv.bind(e); for (size_type i = 0; i < lv.size(); ++i) { auto multiIndex = lv.index(i); if (!compare(U.at(multiIndex), V.at(multiIndex))) return false; } } return true; } // compare discrete functions at quadrature points template bool compare(DiscreteFunction const& u, DiscreteFunction const& v) { 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)); } template void test2(Grid& grid, BasisFactory&& basis) { using namespace Dune::Indices; ProblemStat prob("test", grid, FWD(basis)); prob.initialize(INIT_ALL); auto U0 = *prob.solutionVector(); auto U1 = *prob.solutionVector(); auto U2 = *prob.solutionVector(); auto u0 = valueOf(U0); auto u1 = valueOf(U1); auto u2 = valueOf(U2); // 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); // 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); 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; AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); auto expr2 = [](auto const& x) { return Range{ 1 + 2*x[0] - 3*x[1], 2 + 2*x[0] - 3*x[1]}; }; u0.interpolate_noalias(u2 + expr2); u1.interpolate(u1 + expr2); u2 += expr2; AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); auto expr3 = [](auto const& x) { return Range{ 1 - 0.5*x[0] - 2*x[1], 2 - 0.5*x[0] - 2*x[1]}; }; u0.interpolate_noalias(u2 - expr3); u1.interpolate(u1 - expr3); u2 -= expr3; AMDIS_TEST( compare(U0, U1) ); AMDIS_TEST( compare(U0, U2) ); auto du0 = derivativeOf(u0, tag::gradient{}); auto lf0 = localFunction(u0); auto lf1 = localFunction(u1); auto dlf0 = derivativeOf(lf0, tag::gradient{}); auto dlf1 = derivative(lf1); for (auto const& e : elements(prob.gridView())) { lf0.bind(e); auto geo = e.geometry(); auto local = referenceElement(geo).position(0,0); auto y0 = lf0(local); auto y1 = u0(geo.global(local)); for (auto it1 = y0.begin(), it2 = y1.begin(); it1 != y0.end(); ++it1, ++it2) AMDIS_TEST(compare(*it1, *it2)); // 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(); // should be bound independently to lf0 auto y1_copy = lf0_copy(local); lf0_copy.unbind(); dlf0.bind(e); auto g0 = dlf0(local); auto g1 = du0(geo.global(local)); for (auto it1 = g0.begin(), it2 = g1.begin(); it1 != g0.end(); ++it1, ++it2) AMDIS_TEST(compare(*it1, *it2)); // 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(); } auto V0 = makeDOFVector(prob.globalBasis()); auto v0 = valueOf(V0); v0 << expr1; // test DiscreteFunction std::integral_constant _0; auto tp = makeTreePath(0); auto W0 = *prob.solutionVector(); auto W1 = *prob.solutionVector(); auto W2 = *prob.solutionVector(); auto w0 = valueOf(W0, 0); auto w1 = valueOf(W1, _0); auto w2 = valueOf(W2, tp); // test DiscreteFunction with (pre)treepath argument auto expr = [](auto const& x) { return 1 + x[0] + x[1]; }; auto W3 = *prob.solutionVector(); auto W4 = *prob.solutionVector(); auto W5 = *prob.solutionVector(); auto W7 = *prob.solutionVector(); auto W8 = *prob.solutionVector(); auto w3 = valueOf(W3, 0); auto w4 = valueOf(W4, _0); auto w5 = valueOf(W5, tp); auto w6 = prob.solution(tp); auto& W6 = *prob.solutionVector(); w3 << expr; w4 << expr; w5 << expr; w6 << expr; AMDIS_TEST( compare(W3, W4) ); AMDIS_TEST( compare(W3, W5) ); AMDIS_TEST( compare(W3, W6) ); // test interpolation on subbasis auto w7 = valueOf(W7); auto w8_0 = valueOf(W8, 0); auto w8_1 = valueOf(W8, 1); w7 << expr; w8_0 << expr; w8_1 << expr; 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()); return report_errors(); }