#pragma once #include #include #include #include #include namespace AMDiS { namespace Impl { template auto integrateImpl(GF&& gf, GV const& gridView, QP makeQuad) { auto localFct = localFunction(FWD(gf)); using GridFct = remove_cvref_t; using Range = typename GridFct::Range; Range result(0); for (auto const& element : elements(gridView, Dune::Partitions::interior)) { auto geometry = element.geometry(); localFct.bind(element); auto const& quad = makeQuad(element.type(), localFct); for (auto const qp : quad) result += localFct(qp.position()) * geometry.integrationElement(qp.position()) * qp.weight(); localFct.unbind(); } return gridView.comm().sum(result); } } // end namespace Impl /// \brief Integrate expression with quadrature rule determined by polynomial order of expression /** * **Example:** * ``` * double i1 = integrate(prob.solution(0), prob.gridView()); * double i2 = integrate(X(0) + X(1), prob.gridView()); * double i3 = integrate([](auto const& x) { return x[0] + x[1]; }, prob.gridView()); // ERROR * ``` **/ template auto integrate(Expr&& expr, GridView const& gridView) { auto&& gridFct = makeGridFunction(FWD(expr), gridView); using LocalFct = TYPEOF(localFunction(gridFct)); static_assert(Concepts::Polynomial, "Polynomial degree of expression can not be deduced. You need to provide an explicit value for the quadrature degree or a quadrature rule in `integrate()`."); using Rules = Dune::QuadratureRules; if constexpr(Concepts::Polynomial) return Impl::integrateImpl(FWD(gridFct), gridView, [](auto&& t, auto&& lf) { return Rules::rule(t, order(lf)); }); else return 0.0; } /// \brief Integrate expression with quadrature rule provided /** * **Example:** * ``` * auto quad = Dune::QuadratureRules::rule(Dune::GeometryTypes::triangle, 1); * double i4 = integrate([](auto const& x) { return x[0]; }, prob.gridView(), quad); // OK * ``` **/ template > auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad) { auto&& gridFct = makeGridFunction(FWD(expr), gridView); return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto&&, auto&&) { return quad; }); } /// \brief Integrate expression with quadrature rule determined by provided polynomial `degree` /** * **Example:** * ``` * double i5 = integrate([](auto const& x) { return x[0]; }, prob.gridView(), 1); // OK * ``` **/ template auto integrate(Expr&& expr, GridView const& gridView, int degree, Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre) { using Rules = Dune::QuadratureRules; auto&& gridFct = makeGridFunction(FWD(expr), gridView); return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto const& type, auto&&) { return Rules::rule(type, degree, qt); }); } } // end namespace AMDiS