diff --git a/test/common/CMakeLists.txt b/test/common/CMakeLists.txt index 6670d7e49dc673d5316d8fa09632e60bb5291013..06908ceacc611511f207b98193aea64b87973ba2 100644 --- a/test/common/CMakeLists.txt +++ b/test/common/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(boundingboxtree) +add_subdirectory(functions) add_subdirectory(geometry) add_subdirectory(math) add_subdirectory(parameters) diff --git a/test/common/functions/CMakeLists.txt b/test/common/functions/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..aa9933dd1e1b4cc07f0737375a81111514536e03 --- /dev/null +++ b/test/common/functions/CMakeLists.txt @@ -0,0 +1,3 @@ +dune_add_test(SOURCES test_function_l2norm.cc + CMAKE_GUARD dune-functions_FOUND + LABELS unit) diff --git a/test/common/functions/test_function_l2norm.cc b/test/common/functions/test_function_l2norm.cc new file mode 100644 index 0000000000000000000000000000000000000000..f4911dda842aad4f887ce0be81ca3b5b7862d90c --- /dev/null +++ b/test/common/functions/test_function_l2norm.cc @@ -0,0 +1,80 @@ +#include <config.h> + +#include <dune/common/fvector.hh> +#include <dune/istl/bvector.hh> + +#include <dune/functions/functionspacebases/interpolate.hh> +#include <dune/functions/functionspacebases/lagrangebasis.hh> +#include <dune/functions/gridfunctions/analyticgridviewfunction.hh> +#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> + +#include <dumux/common/parameters.hh> +#include <dumux/common/integrate.hh> +#include <dumux/multidomain/glue.hh> +#include <dumux/discretization/projection/projector.hh> +#include <dumux/discretization/fem/fegridgeometry.hh> +#include <dumux/io/grid/gridmanager_yasp.hh> + +int main (int argc, char *argv[]) try +{ + // maybe initialize mpi + Dune::MPIHelper::instance(argc, argv); + + // initialize parameters + Dumux::Parameters::init([](auto& p){ + p["Grid.UpperRight"] = "1 1"; + p["Coarse.Grid.Cells"] = "10 10"; + p["Fine.Grid.Cells"] = "33 33"; + }); + + // initialize a grid + using Grid = Dune::YaspGrid<2>; + Dumux::GridManager<Grid> gridManagerCoarse, gridManagerFine; + gridManagerCoarse.init("Coarse"); + gridManagerFine.init("Fine"); + const auto gvFine = gridManagerFine.grid().leafGridView(); + const auto gvCoarse = gridManagerCoarse.grid().leafGridView(); + + // a quadratic function on the grid + auto f = [](const auto& pos){ return pos.two_norm2(); }; + + // make discrete functions + using R = Dune::FieldVector<double, 1>; + using SolutionVector = Dune::BlockVector<R>; + using GridView = Grid::LeafGridView; + using namespace Dune::Functions; + using P2 = LagrangeBasis<GridView, 2>; + auto basisCoarse = std::make_shared<P2>(gvCoarse), basisFine = std::make_shared<P2>(gvFine); + + SolutionVector solCoarse, solFine; + interpolate(*basisCoarse, solCoarse, f); + interpolate(*basisFine, solFine, f); + auto gfCoarse = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse); + auto gfFine = makeDiscreteGlobalBasisFunction<R>(*basisFine, solFine); + + // project fine discrete function onto coarse grid + Dumux::FEGridGeometry<P2> ggCoarse(basisCoarse), ggFine(basisFine); + const auto projector = Dumux::makeProjector(*basisFine, *basisCoarse, makeGlue(ggFine, ggCoarse)); + SolutionVector solCoarse2; + projector.project(solFine, solCoarse2); + auto gfCoarse2 = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse); + + // check that these functions are identical and that they exactly represent the analytical function + auto l2norm = Dumux::integrateL2Error(gvCoarse, gfCoarse, gfCoarse2, 3); + if (l2norm > 1e-14) + DUNE_THROW(Dune::MathError, "L2 norm of projection and interpolation is non-zero!"); + + l2norm = Dumux::integrateL2Error(gvCoarse, makeAnalyticGridViewFunction(f, gvCoarse), gfCoarse2, 3); + if (l2norm > 1e-14) + DUNE_THROW(Dune::MathError, "L2 norm of analytical solution and interpolation is non-zero!"); + + std::cout << "All tests passed!" << std::endl; + return 0; +} +// ////////////////////////////////// +// Error handler +// ///////////////////////////////// +catch (const Dune::Exception& e) { + std::cout << e << std::endl; + return 1; +}