From e08eab8a36bc2a56548c181109bbfe4234cdc019 Mon Sep 17 00:00:00 2001 From: Felix Schindler <felix.schindler@wwu.de> Date: Sun, 10 Apr 2016 13:38:28 +0200 Subject: [PATCH] [linearelliptic.estimators.swipdg...] refactor to match current state --- .../estimators/swipdg-fluxreconstruction.hh | 1219 +++++++---------- 1 file changed, 467 insertions(+), 752 deletions(-) diff --git a/dune/gdt/test/linearelliptic/estimators/swipdg-fluxreconstruction.hh b/dune/gdt/test/linearelliptic/estimators/swipdg-fluxreconstruction.hh index 97949bfc5..615e9c16d 100644 --- a/dune/gdt/test/linearelliptic/estimators/swipdg-fluxreconstruction.hh +++ b/dune/gdt/test/linearelliptic/estimators/swipdg-fluxreconstruction.hh @@ -1,966 +1,681 @@ -// This file is part of the dune-hdd project: -// http://users.dune-project.org/projects/dune-hdd +// This file is part of the dune-gdt project: +// http://users.dune-project.org/projects/dune-gdt // Copyright holders: Felix Schindler // License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) -#ifndef DUNE_HDD_LINEARELLIPTIC_ESTIMATORS_SWIPDG_HH -#define DUNE_HDD_LINEARELLIPTIC_ESTIMATORS_SWIPDG_HH +#ifndef DUNE_GDT_TEST_LINEARELLIPTIC_ESTIMATORS_SWIPDG_FLUXRECONSTRUCTION_HH +#define DUNE_GDT_TEST_LINEARELLIPTIC_ESTIMATORS_SWIPDG_FLUXRECONSTRUCTION_HH -#include <map> -#include <memory> -#include <string> -#include <vector> - -#include <boost/numeric/conversion/cast.hpp> - -#if HAVE_ALUGRID -#include <dune/grid/alugrid.hh> -#endif - -#include <dune/stuff/common/memory.hh> -#include <dune/stuff/common/tmp-storage.hh> #include <dune/stuff/grid/walker.hh> #include <dune/stuff/grid/walker/functors.hh> #include <dune/stuff/playground/functions/ESV2007.hh> - -#include <dune/pymor/common/exceptions.hh> -#include <dune/pymor/parameters/base.hh> +#include <dune/stuff/la/container.hh> #include <dune/gdt/discretefunction/default.hh> #include <dune/gdt/localevaluation/elliptic.hh> #include <dune/gdt/localevaluation/product.hh> -#include <dune/gdt/localoperator/codim0.hh> +#include <dune/gdt/localoperator/integrals.hh> #include <dune/gdt/operators/oswaldinterpolation.hh> -#include <dune/gdt/operators/projections.hh> +#include <dune/gdt/projections.hh> #include <dune/gdt/playground/localevaluation/ESV2007.hh> #include <dune/gdt/playground/operators/fluxreconstruction.hh> -#include <dune/gdt/playground/spaces/finitevolume/default.hh> +#include <dune/gdt/spaces/fv/default.hh> +#include <dune/gdt/spaces/rt/pdelab.hh> namespace Dune { -namespace HDD { +namespace GDT { namespace LinearElliptic { -namespace Estimators { -namespace internal { -namespace SWIPDG { +namespace SwipdgFluxreconstrutionEstimators { static const size_t over_integrate = 2; -class LocalNonconformityESV2007Base +static std::string local_nonconformity_ESV2007_id() { -public: - static std::string id() - { - return "eta_NC_ESV2007"; - } -}; + return "eta_NC_ESV2007"; +} -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class LocalNonconformityESV2007 : public LocalNonconformityESV2007Base +static std::string local_residual_ESV2007_id() { -public: - static const bool available = false; -}; + return "eta_R_ESV2007"; +} -#if HAVE_ALUGRID -/** - * \brief computes the local nonconformity estimator as defined in ESV2007 - */ -template <class SpaceType, class VectorType, class ProblemType> -class LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> - : public LocalNonconformityESV2007Base, public Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> +static std::string local_diffusive_flux_ESV2007_id() { - typedef LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> ThisType; - typedef Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> FunctorBaseType; - -public: - static const bool available = true; - - typedef std::map<std::string, Pymor::Parameter> ParametersMapType; - - typedef typename FunctorBaseType::GridViewType GridViewType; - typedef typename FunctorBaseType::EntityType EntityType; + return "eta_DF_ESV2007"; +} - typedef typename ProblemType::RangeFieldType RangeFieldType; -private: - typedef GDT::ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; - typedef GDT::DiscreteFunction<SpaceType, VectorType> DiscreteFunctionType; - typedef typename ConstDiscreteFunctionType::DifferenceType DifferenceType; - - typedef typename ProblemType::DiffusionFactorType::NonparametricType DiffusionFactorType; - typedef typename ProblemType::DiffusionTensorType::NonparametricType DiffusionTensorType; - - typedef GDT::LocalOperator::Codim0Integral<GDT::LocalEvaluation::Elliptic<DiffusionFactorType, DiffusionTensorType>> - LocalOperatorType; - typedef Stuff::Common::TmpMatricesStorage<RangeFieldType> TmpStorageProviderType; - - static const ProblemType& assert_problem(const ProblemType& problem, const Pymor::Parameter& mu_bar) - { - if (mu_bar.type() != problem.parameter_type()) - DUNE_THROW(Pymor::Exceptions::wrong_parameter_type, - "Given mu_bar is of type " << mu_bar.type() << " and should be of type " << problem.parameter_type() - << "!"); - if (problem.diffusion_tensor()->parametric()) - DUNE_THROW(NotImplemented, "Not implemented for parametric diffusion_tensor!"); - return problem; - } // ... assert_problem(...) - -public: - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const ParametersMapType parameters = ParametersMapType()) - { - if (problem.diffusion_factor()->parametric() && parameters.find("mu_bar") == parameters.end()) - DUNE_THROW(Stuff::Exceptions::wrong_input_given, "Given parameters are missing 'mu_bar'!"); - const Pymor::Parameter mu_bar = problem.parametric() ? parameters.at("mu_bar") : Pymor::Parameter(); - ThisType estimator(space, vector, problem, mu_bar); - Stuff::Grid::Walker<GridViewType> grid_walker(space.grid_view()); - grid_walker.add(estimator); - grid_walker.walk(); - return std::sqrt(estimator.result_); - } // ... estimate(...) - - LocalNonconformityESV2007(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const Pymor::Parameter mu_bar = Pymor::Parameter()) - : space_(space) - , vector_(vector) - , problem_(assert_problem(problem, mu_bar)) - , problem_mu_bar_(problem_.with_mu(mu_bar)) - , discrete_solution_(space_, vector_) - , oswald_interpolation_(space_) - , difference_(Stuff::Common::make_unique<DifferenceType>(discrete_solution_ - oswald_interpolation_)) - , local_operator_(over_integrate, *problem_mu_bar_->diffusion_factor()->affine_part(), - *problem_.diffusion_tensor()->affine_part()) - , tmp_local_matrices_({1, local_operator_.numTmpObjectsRequired()}, 1, 1) - , prepared_(false) - , result_(0.0) - { - } - - virtual void prepare() - { - if (!prepared_) { - const GDT::Operators::OswaldInterpolation<GridViewType> oswald_interpolation_operator(space_.grid_view()); - oswald_interpolation_operator.apply(discrete_solution_, oswald_interpolation_); - result_ = 0.0; - prepared_ = true; - } - } // ... prepare(...) - - RangeFieldType compute_locally(const EntityType& entity) - { - const auto local_difference = difference_->local_function(entity); - local_operator_.apply( - *local_difference, *local_difference, tmp_local_matrices_.matrices()[0][0], tmp_local_matrices_.matrices()[1]); - assert(tmp_local_matrices_.matrices()[0][0].rows() >= 1); - assert(tmp_local_matrices_.matrices()[0][0].cols() >= 1); - return tmp_local_matrices_.matrices()[0][0][0][0]; - } // ... compute_locally(...) - - virtual void apply_local(const EntityType& entity) - { - result_ += compute_locally(entity); - } - -private: - const SpaceType& space_; - const VectorType& vector_; - const ProblemType& problem_; - const std::shared_ptr<const typename ProblemType::NonparametricType> problem_mu_bar_; - const ConstDiscreteFunctionType discrete_solution_; - DiscreteFunctionType oswald_interpolation_; - std::unique_ptr<const DifferenceType> difference_; - const LocalOperatorType local_operator_; - TmpStorageProviderType tmp_local_matrices_; - bool prepared_; - -public: - RangeFieldType result_; -}; // class LocalNonconformityESV2007< ..., ALUGrid< 2, 2, simplex, conforming >, ... > - -#endif // HAVE_ALUGRID - - -class LocalResidualESV2007Base +static std::string ESV2007_id() { -public: - static std::string id() - { - return "eta_R_ESV2007"; - } -}; + return "eta_ESV2007"; +} -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class LocalResidualESV2007 : public LocalResidualESV2007Base +static std::string ESV2007_alternative_summation_id() { -public: - static const bool available = false; -}; + return ESV2007_id() + "_alternative_summation"; +} -#if HAVE_ALUGRID /** - * \brief computes the local residual estimator as defined in ESV2007 + * \brief computes the local nonconformity estimator as defined in ESV2007 */ -template <class SpaceType, class VectorType, class ProblemType> -class LocalResidualESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> - : public LocalResidualESV2007Base, public Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> +template <class SpaceType, class VectorType, class DiffusionFactorType, class DiffusionTensorType, + class GridViewType = typename SpaceType::GridViewType> +class LocalNonconformityESV2007 + : public Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename SpaceType::RangeFieldType> { - typedef LocalResidualESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> ThisType; - typedef Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> FunctorBaseType; + typedef Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename SpaceType::RangeFieldType> BaseType; + typedef LocalNonconformityESV2007<SpaceType, VectorType, DiffusionFactorType, DiffusionTensorType, GridViewType> + ThisType; + typedef ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; + typedef DiscreteFunction<SpaceType, VectorType> DiscreteFunctionType; + typedef typename ConstDiscreteFunctionType::DifferenceType DifferenceType; + typedef LocalVolumeIntegralOperator<LocalEvaluation::Elliptic<DiffusionFactorType, DiffusionTensorType>> + LocalOperatorType; public: - static const bool available = true; - - typedef typename FunctorBaseType::GridViewType GridViewType; - typedef typename FunctorBaseType::EntityType EntityType; - - typedef typename ProblemType::RangeFieldType RangeFieldType; + using typename BaseType::EntityType; + using typename BaseType::ReturnType; -private: - typedef GDT::Spaces::FiniteVolume::Default<GridViewType, RangeFieldType, 1, 1> P0SpaceType; - typedef GDT::DiscreteFunction<P0SpaceType, VectorType> DiscreteFunctionType; - typedef typename DiscreteFunctionType::DifferenceType DifferenceType; - - typedef typename ProblemType::DiffusionFactorType::NonparametricType DiffusionFactorType; - typedef typename ProblemType::DiffusionTensorType::NonparametricType DiffusionTensorType; - - typedef typename Stuff::Functions::ESV2007::Cutoff<DiffusionFactorType, DiffusionTensorType> CutoffFunctionType; - typedef GDT::LocalOperator::Codim0Integral<GDT::LocalEvaluation::Product<CutoffFunctionType>> LocalOperatorType; - typedef DSC::TmpMatricesStorage<RangeFieldType> TmpStorageProviderType; + static std::string id() + { + return local_nonconformity_ESV2007_id(); + } - static const ProblemType& assert_problem(const ProblemType& problem) + static ReturnType estimate(const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) { - if (problem.parametric()) - DUNE_THROW(NotImplemented, "Not implemented yet for parametric problems!"); - assert(problem.diffusion_factor()->has_affine_part()); - assert(problem.diffusion_tensor()->has_affine_part()); - assert(problem.force()->has_affine_part()); - return problem; - } // ... assert_problem(...) + return estimate(space.grid_view(), space, vector, diffusion_factor, diffusion_tensor, over_int); + } -public: - static RangeFieldType estimate(const SpaceType& space, const VectorType& /*vector*/, const ProblemType& problem) + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) { - ThisType estimator(space, problem); - Stuff::Grid::Walker<GridViewType> grid_walker(space.grid_view()); + ThisType estimator(grid_view, space, vector, diffusion_factor, diffusion_tensor, over_int); + Stuff::Grid::Walker<GridViewType> grid_walker(grid_view); grid_walker.add(estimator); grid_walker.walk(); - return std::sqrt(estimator.result_); + return std::sqrt(estimator.result()); } // ... estimate(...) - LocalResidualESV2007(const SpaceType& space, const ProblemType& problem) - : space_(space) - , problem_(assert_problem(problem)) - , p0_space_(space_.grid_view()) - , p0_force_(p0_space_) - , difference_(Stuff::Common::make_unique<DifferenceType>(*problem_.force()->affine_part() - p0_force_)) - , cutoff_function_(*problem_.diffusion_factor()->affine_part(), *problem_.diffusion_tensor()->affine_part()) - , local_operator_(over_integrate, cutoff_function_) - , tmp_local_matrices_({1, local_operator_.numTmpObjectsRequired()}, 1, 1) + LocalNonconformityESV2007(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) + : grid_view_(grid_view) + , discrete_solution_(space, vector) + , oswald_interpolation_(space) + , difference_(discrete_solution_ - oswald_interpolation_) + , local_operator_(over_int, diffusion_factor, diffusion_tensor) , prepared_(false) , result_(0.0) { } - virtual void prepare() + virtual void prepare() override final { - if (!prepared_) { - const GDT::Operators::Projection<GridViewType> projection_operator(space_.grid_view(), over_integrate); - projection_operator.apply(*problem_.force()->affine_part(), p0_force_); - result_ = 0.0; - prepared_ = true; - } + if (prepared_) + return; + const Operators::OswaldInterpolation<GridViewType> oswald_interpolation_operator(grid_view_); + oswald_interpolation_operator.apply(discrete_solution_, oswald_interpolation_); + result_ = 0.0; + prepared_ = true; } // ... prepare(...) - RangeFieldType compute_locally(const EntityType& entity) + virtual ReturnType compute_locally(const EntityType& entity) override final { - const auto local_difference = difference_->local_function(entity); - local_operator_.apply( - *local_difference, *local_difference, tmp_local_matrices_.matrices()[0][0], tmp_local_matrices_.matrices()[1]); - assert(tmp_local_matrices_.matrices()[0][0].rows() >= 1); - assert(tmp_local_matrices_.matrices()[0][0].cols() >= 1); - return tmp_local_matrices_.matrices()[0][0][0][0]; + const auto local_difference = difference_.local_function(entity); + const auto result = local_operator_.apply2(*local_difference, *local_difference); + assert(result.rows() >= 1); + assert(result.cols() >= 1); + return result[0][0]; } // ... compute_locally(...) - virtual void apply_local(const EntityType& entity) + virtual void apply_local(const EntityType& entity) override final { result_ += compute_locally(entity); } -private: - const SpaceType& space_; - const ProblemType& problem_; - const P0SpaceType p0_space_; - DiscreteFunctionType p0_force_; - std::unique_ptr<const DifferenceType> difference_; - const CutoffFunctionType cutoff_function_; - const LocalOperatorType local_operator_; - TmpStorageProviderType tmp_local_matrices_; - bool prepared_; - -public: - RangeFieldType result_; -}; // class LocalResidualESV2007< ..., ALUGrid< 2, 2, simplex, conforming >, ... > - -#endif // HAVE_ALUGRID - - -class LocalResidualESV2007StarBase -{ -public: - static std::string id() + virtual ReturnType result() const { - return "eta_R_ESV2007_*"; + return result_; } -}; +private: + const GridViewType& grid_view_; + const ConstDiscreteFunctionType discrete_solution_; + DiscreteFunctionType oswald_interpolation_; + const DifferenceType difference_; + const LocalOperatorType local_operator_; + bool prepared_; + ReturnType result_; +}; // class LocalNonconformityESV2007 -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class LocalResidualESV2007Star : public LocalResidualESV2007StarBase -{ -public: - static const bool available = false; -}; - -#if HAVE_ALUGRID /** * \brief computes the local residual estimator as defined in ESV2007 */ -template <class SpaceType, class VectorType, class ProblemType> -class LocalResidualESV2007Star<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> - : public LocalResidualESV2007StarBase, public Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> +template <class SpaceType, class VectorType, class ForceType, class DiffusionFactorType, class DiffusionTensorType, + class GridViewType = typename SpaceType::GridViewType> +class LocalResidualESV2007 : public Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename ForceType::RangeFieldType> { - typedef LocalResidualESV2007Star<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> ThisType; - typedef Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> FunctorBaseType; - -public: - static const bool available = true; - - typedef std::map<std::string, Pymor::Parameter> ParametersMapType; - - typedef typename FunctorBaseType::GridViewType GridViewType; - typedef typename FunctorBaseType::EntityType EntityType; - - typedef typename ProblemType::RangeFieldType RangeFieldType; - - static const unsigned int dimDomain = SpaceType::dimDomain; - -private: - typedef GDT::ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; - typedef GDT::Spaces::RaviartThomas::PdelabBased<GridViewType, 0, RangeFieldType, dimDomain> RTN0SpaceType; - typedef GDT::DiscreteFunction<RTN0SpaceType, VectorType> RTN0DiscreteFunctionType; - typedef typename RTN0DiscreteFunctionType::DivergenceType DivergenceType; + typedef Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename ForceType::RangeFieldType> BaseType; + typedef LocalResidualESV2007<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, GridViewType> + ThisType; + typedef typename ForceType::RangeFieldType RangeFieldType; + typedef ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; + typedef Spaces::RT::PdelabBased<GridViewType, 0, RangeFieldType, SpaceType::dimDomain> RTN0SpaceType; + typedef DiscreteFunction<RTN0SpaceType, VectorType> DiffusiveFluxType; + typedef typename DiffusiveFluxType::DivergenceType DivergenceType; typedef typename DivergenceType::DifferenceType DifferenceType; + typedef typename Stuff::Functions::ESV2007::Cutoff<DiffusionFactorType, DiffusionTensorType> CutoffFunctionType; + typedef LocalVolumeIntegralOperator<LocalEvaluation::Product<CutoffFunctionType>> LocalOperatorType; - typedef typename ProblemType::DiffusionFactorType::NonparametricType DiffusionFactorType; - typedef typename ProblemType::DiffusionTensorType::NonparametricType DiffusionTensorType; +public: + using typename BaseType::EntityType; + using typename BaseType::ReturnType; - typedef typename Stuff::Functions::ESV2007::Cutoff<DiffusionFactorType, DiffusionTensorType> CutoffFunctionType; - typedef GDT::LocalOperator::Codim0Integral<GDT::LocalEvaluation::Product<CutoffFunctionType>> LocalOperatorType; - typedef DSC::TmpMatricesStorage<RangeFieldType> TmpStorageProviderType; + static std::string id() + { + return local_residual_ESV2007_id(); + } - static const ProblemType& assert_problem(const ProblemType& problem, const Pymor::Parameter& mu) + static ReturnType estimate(const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) { - if (mu.type() != problem.parameter_type()) - DUNE_THROW(Pymor::Exceptions::wrong_parameter_type, - "Given mu is of type " << mu.type() << " and should be of type " << problem.parameter_type() << "!"); - if (problem.diffusion_tensor()->parametric()) - DUNE_THROW(NotImplemented, "Not implemented for parametric diffusion_tensor!"); - if (problem.force()->parametric()) - DUNE_THROW(NotImplemented, "Not implemented for parametric force!"); - return problem; - } // ... assert_problem(...) + return estimate(space.grid_view(), space, vector, force, diffusion_factor, diffusion_tensor, over_int); + } -public: - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const ParametersMapType parameters = ParametersMapType()) - { - if (problem.diffusion_factor()->parametric() && parameters.find("mu") == parameters.end()) - DUNE_THROW(Stuff::Exceptions::wrong_input_given, "Given parameters are missing 'mu'!"); - const Pymor::Parameter mu = problem.parametric() ? parameters.at("mu") : Pymor::Parameter(); - ThisType estimator(space, vector, problem, mu); - Stuff::Grid::Walker<GridViewType> grid_walker(space.grid_view()); + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + ThisType estimator(grid_view, space, vector, force, diffusion_factor, diffusion_tensor, over_int); + Stuff::Grid::Walker<GridViewType> grid_walker(grid_view); grid_walker.add(estimator); grid_walker.walk(); - return std::sqrt(estimator.result_); + return std::sqrt(estimator.result()); } // ... estimate(...) - LocalResidualESV2007Star(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const Pymor::Parameter mu = Pymor::Parameter()) - : space_(space) - , vector_(vector) - , problem_(assert_problem(problem, mu)) - , problem_mu_(problem_.with_mu(mu)) - , discrete_solution_(space_, vector_) - , rtn0_space_(space_.grid_view()) + LocalResidualESV2007(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + : grid_view_(grid_view) + , diffusion_factor_(diffusion_factor) + , diffusion_tensor_(diffusion_tensor) + , over_integrate_(over_int) + , discrete_solution_(space, vector) + , rtn0_space_(grid_view_) , diffusive_flux_(rtn0_space_) , divergence_(diffusive_flux_.divergence()) - , difference_(*problem_.force()->affine_part() - divergence_) - , cutoff_function_(*problem_.diffusion_factor()->affine_part(), *problem_.diffusion_tensor()->affine_part()) - , local_operator_(over_integrate, cutoff_function_) - , tmp_local_matrices_({1, local_operator_.numTmpObjectsRequired()}, 1, 1) + , difference_(force - divergence_) + , cutoff_function_(diffusion_factor_, diffusion_tensor_) + , local_operator_(over_integrate_, cutoff_function_) , prepared_(false) , result_(0.0) { } - virtual ~LocalResidualESV2007Star() = default; - - virtual void prepare() + virtual void prepare() override final { - if (!prepared_) { - const GDT::Operators::DiffusiveFluxReconstruction<GridViewType, DiffusionFactorType, DiffusionTensorType> - diffusive_flux_reconstruction(space_.grid_view(), - *problem_mu_->diffusion_factor()->affine_part(), - *problem_.diffusion_tensor()->affine_part(), - over_integrate); - diffusive_flux_reconstruction.apply(discrete_solution_, diffusive_flux_); - result_ = 0.0; - prepared_ = true; - } + if (prepared_) + return; + const Operators::DiffusiveFluxReconstruction<GridViewType, DiffusionFactorType, DiffusionTensorType> + diffusive_flux_reconstruction(grid_view_, diffusion_factor_, diffusion_tensor_, over_integrate_); + diffusive_flux_reconstruction.apply(discrete_solution_, diffusive_flux_); + result_ = 0.0; + prepared_ = true; } // ... prepare(...) - RangeFieldType compute_locally(const EntityType& entity) + virtual ReturnType compute_locally(const EntityType& entity) override final { const auto local_difference = difference_.local_function(entity); - local_operator_.apply( - *local_difference, *local_difference, tmp_local_matrices_.matrices()[0][0], tmp_local_matrices_.matrices()[1]); - assert(tmp_local_matrices_.matrices()[0][0].rows() >= 1); - assert(tmp_local_matrices_.matrices()[0][0].cols() >= 1); - return tmp_local_matrices_.matrices()[0][0][0][0]; + auto result = local_operator_.apply2(*local_difference, *local_difference); + assert(result.rows() >= 1); + assert(result.cols() >= 1); + return result[0][0]; } // ... compute_locally(...) - virtual void apply_local(const EntityType& entity) + virtual void apply_local(const EntityType& entity) override final { result_ += compute_locally(entity); } + virtual ReturnType result() const override final + { + return result_; + } + private: - const SpaceType& space_; - const VectorType& vector_; - const ProblemType& problem_; - const std::shared_ptr<typename ProblemType::NonparametricType> problem_mu_; + const GridViewType grid_view_; + const DiffusionFactorType& diffusion_factor_; + const DiffusionTensorType& diffusion_tensor_; + const size_t over_integrate_; const ConstDiscreteFunctionType discrete_solution_; const RTN0SpaceType rtn0_space_; - RTN0DiscreteFunctionType diffusive_flux_; + DiffusiveFluxType diffusive_flux_; const DivergenceType divergence_; const DifferenceType difference_; const CutoffFunctionType cutoff_function_; const LocalOperatorType local_operator_; - TmpStorageProviderType tmp_local_matrices_; bool prepared_; - -public: RangeFieldType result_; -}; // class LocalResidualESV2007Star< ..., ALUGrid< 2, 2, simplex, conforming >, ... > +}; // class LocalResidualESV2007 -#endif // HAVE_ALUGRID - - -class LocalDiffusiveFluxESV2007Base -{ -public: - static std::string id() - { - return "eta_DF_ESV2007"; - } -}; - - -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class LocalDiffusiveFluxESV2007 : public LocalDiffusiveFluxESV2007Base -{ -public: - static const bool available = false; -}; - -#if HAVE_ALUGRID /** * \brief computes the local diffusive flux estimator as defined in ESV2007 */ -template <class SpaceType, class VectorType, class ProblemType> -class LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> - : public LocalDiffusiveFluxESV2007Base, public Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> +template <class SpaceType, class VectorType, class DiffusionFactorType, class DiffusionTensorType, + class GridViewType = typename SpaceType::GridViewType> +class LocalDiffusiveFluxESV2007 + : public Stuff::Grid::Codim0ReturnFunctor<typename SpaceType::GridViewType, typename SpaceType::RangeFieldType> { - typedef LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> ThisType; - typedef Stuff::Grid::Functor::Codim0<typename SpaceType::GridViewType> FunctorBaseType; + typedef LocalDiffusiveFluxESV2007<SpaceType, VectorType, DiffusionFactorType, DiffusionTensorType, GridViewType> + ThisType; + typedef Stuff::Grid::Codim0ReturnFunctor<typename SpaceType::GridViewType, typename SpaceType::RangeFieldType> + BaseType; + typedef typename SpaceType::RangeFieldType RangeFieldType; + typedef ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; + typedef Spaces::RT::PdelabBased<GridViewType, 0, RangeFieldType, SpaceType::dimDomain> RTN0SpaceType; + typedef DiscreteFunction<RTN0SpaceType, VectorType> RTN0DiscreteFunctionType; + typedef LocalVolumeIntegralOperator<LocalEvaluation::ESV2007::DiffusiveFluxEstimate<DiffusionFactorType, + RTN0DiscreteFunctionType, + DiffusionTensorType>> + LocalOperatorType; public: - static const bool available = true; - - typedef std::map<std::string, Pymor::Parameter> ParametersMapType; + using typename BaseType::EntityType; + using typename BaseType::ReturnType; - typedef typename FunctorBaseType::GridViewType GridViewType; - typedef typename FunctorBaseType::EntityType EntityType; + static std::string id() + { + return local_diffusive_flux_ESV2007_id(); + } - typedef typename ProblemType::RangeFieldType RangeFieldType; + static ReturnType estimate(const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), space, vector, diffusion_factor, diffusion_tensor, over_int); + } - static const unsigned int dimDomain = SpaceType::dimDomain; + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) + { + return estimate(grid_view, space, vector, diffusion_factor, diffusion_factor, diffusion_tensor, over_int); + } -private: - typedef GDT::ConstDiscreteFunction<SpaceType, VectorType> ConstDiscreteFunctionType; - typedef GDT::Spaces::RaviartThomas::PdelabBased<GridViewType, 0, RangeFieldType, dimDomain> RTN0SpaceType; - typedef GDT::DiscreteFunction<RTN0SpaceType, VectorType> RTN0DiscreteFunctionType; - - typedef typename ProblemType::DiffusionFactorType::NonparametricType DiffusionFactorType; - typedef typename ProblemType::DiffusionTensorType::NonparametricType DiffusionTensorType; - - typedef GDT::LocalOperator:: - Codim0Integral<GDT::LocalEvaluation::ESV2007::DiffusiveFluxEstimate<DiffusionFactorType, RTN0DiscreteFunctionType, - DiffusionTensorType>> LocalOperatorType; - typedef DSC::TmpMatricesStorage<RangeFieldType> TmpStorageProviderType; - - static const ProblemType& assert_problem(const ProblemType& problem, const Pymor::Parameter& mu, - const Pymor::Parameter& mu_hat) - { - if (mu.type() != problem.parameter_type()) - DUNE_THROW(Pymor::Exceptions::wrong_parameter_type, - "Given mu is of type " << mu.type() << " and should be of type " << problem.parameter_type() << "!"); - if (mu_hat.type() != problem.parameter_type()) - DUNE_THROW(Pymor::Exceptions::wrong_parameter_type, - "Given mu_hat is of type " << mu_hat.type() << " and should be of type " << problem.parameter_type() - << "!"); - if (problem.diffusion_tensor()->parametric()) - DUNE_THROW(NotImplemented, "Not implemented for parametric diffusion_tensor!"); - return problem; - } // ... assert_problem(...) + static ReturnType estimate(const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), + space, + vector, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); + } -public: - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const ParametersMapType parameters = ParametersMapType()) - { - if (problem.diffusion_factor()->parametric() && parameters.find("mu") == parameters.end()) - DUNE_THROW(Stuff::Exceptions::wrong_input_given, "Given parameters are missing 'mu'!"); - if (problem.diffusion_factor()->parametric() && parameters.find("mu_hat") == parameters.end()) - DUNE_THROW(Stuff::Exceptions::wrong_input_given, "Given parameters are missing 'mu_hat'!"); - const Pymor::Parameter mu = problem.parametric() ? parameters.at("mu") : Pymor::Parameter(); - const Pymor::Parameter mu_hat = problem.parametric() ? parameters.at("mu_hat") : Pymor::Parameter(); - ThisType estimator(space, vector, problem, mu, mu_hat); - Stuff::Grid::Walker<GridViewType> grid_walker(space.grid_view()); + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + ThisType estimator( + grid_view, space, vector, diffusion_factor_norm, diffusion_factor_reconstruction, diffusion_tensor, over_int); + Stuff::Grid::Walker<GridViewType> grid_walker(grid_view); grid_walker.add(estimator); grid_walker.walk(); - return std::sqrt(estimator.result_); + return std::sqrt(estimator.result()); } // ... estimate(...) - LocalDiffusiveFluxESV2007(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const Pymor::Parameter mu = Pymor::Parameter(), - const Pymor::Parameter mu_hat = Pymor::Parameter()) - : space_(space) - , vector_(vector) - , problem_(assert_problem(problem, mu, mu_hat)) - , problem_mu_(problem_.with_mu(mu)) - , problem_mu_hat_(problem.with_mu(mu_hat)) - , discrete_solution_(space_, vector_) - , rtn0_space_(space.grid_view()) + LocalDiffusiveFluxESV2007(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + : grid_view_(grid_view) + , diffusion_factor_reconstruction_(diffusion_factor_reconstruction) + , diffusion_tensor_(diffusion_tensor) + , discrete_solution_(space, vector) + , rtn0_space_(grid_view_) , diffusive_flux_(rtn0_space_) - , local_operator_(over_integrate, *problem_mu_hat_->diffusion_factor()->affine_part(), - *problem_.diffusion_tensor()->affine_part(), diffusive_flux_) - , tmp_local_matrices_({1, local_operator_.numTmpObjectsRequired()}, 1, 1) + , over_int_(over_int) + , local_operator_(over_int_, diffusion_factor_norm, diffusion_tensor_, diffusive_flux_) , prepared_(false) , result_(0.0) { } - virtual void prepare() + virtual void prepare() override final { - if (!prepared_) { - const GDT::Operators::DiffusiveFluxReconstruction<GridViewType, DiffusionFactorType, DiffusionTensorType> - diffusive_flux_reconstruction(space_.grid_view(), - *problem_mu_->diffusion_factor()->affine_part(), - *problem_.diffusion_tensor()->affine_part(), - over_integrate); - diffusive_flux_reconstruction.apply(discrete_solution_, diffusive_flux_); - result_ = 0.0; - prepared_ = true; - } + if (prepared_) + return; + const Operators::DiffusiveFluxReconstruction<GridViewType, DiffusionFactorType, DiffusionTensorType> + diffusive_flux_reconstruction(grid_view_, diffusion_factor_reconstruction_, diffusion_tensor_, over_int_); + diffusive_flux_reconstruction.apply(discrete_solution_, diffusive_flux_); + result_ = 0.0; + prepared_ = true; } // ... prepare(...) - RangeFieldType compute_locally(const EntityType& entity) + virtual ReturnType compute_locally(const EntityType& entity) override final { const auto local_discrete_solution = discrete_solution_.local_function(entity); - local_operator_.apply(*local_discrete_solution, - *local_discrete_solution, - tmp_local_matrices_.matrices()[0][0], - tmp_local_matrices_.matrices()[1]); - assert(tmp_local_matrices_.matrices()[0][0].rows() >= 1); - assert(tmp_local_matrices_.matrices()[0][0].cols() >= 1); - return tmp_local_matrices_.matrices()[0][0][0][0]; + auto result = local_operator_.apply2(*local_discrete_solution, *local_discrete_solution); + assert(result.rows() >= 1); + assert(result.cols() >= 1); + return result[0][0]; } // ... compute_locally(...) - virtual void apply_local(const EntityType& entity) + virtual void apply_local(const EntityType& entity) override final { result_ += compute_locally(entity); } + virtual ReturnType result() const override final + { + return result_; + } + private: - const SpaceType& space_; - const VectorType& vector_; - const ProblemType& problem_; - const std::shared_ptr<const typename ProblemType::NonparametricType> problem_mu_; - const std::shared_ptr<const typename ProblemType::NonparametricType> problem_mu_hat_; + const GridViewType& grid_view_; + const DiffusionFactorType& diffusion_factor_reconstruction_; + const DiffusionTensorType& diffusion_tensor_; const ConstDiscreteFunctionType discrete_solution_; const RTN0SpaceType rtn0_space_; RTN0DiscreteFunctionType diffusive_flux_; + const size_t over_int_; const LocalOperatorType local_operator_; - TmpStorageProviderType tmp_local_matrices_; bool prepared_; - -public: RangeFieldType result_; -}; // class LocalDiffusiveFluxESV2007< ..., ALUGrid< 2, 2, simplex, conforming >, ... > - -#endif // HAVE_ALUGRID +}; // class LocalDiffusiveFluxESV2007 -class ESV2007Base +template <class SpaceType, class VectorType, class ForceType, class DiffusionFactorType, class DiffusionTensorType, + class GridViewType = typename SpaceType::GridViewType> +class ESV2007 : public Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename SpaceType::RangeFieldType> { -public: - static std::string id() - { - return "eta_ESV2007"; - } -}; - - -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class ESV2007 : public ESV2007Base -{ -public: - static const bool available = false; -}; + typedef Stuff::Grid::Codim0ReturnFunctor<GridViewType, typename SpaceType::RangeFieldType> BaseType; + typedef ESV2007<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, GridViewType> ThisType; - -#if HAVE_ALUGRID - -template <class SpaceType, class VectorType, class ProblemType> -class ESV2007<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> : public ESV2007Base -{ - typedef ALUGrid<2, 2, simplex, conforming> GridType; + typedef LocalNonconformityESV2007<SpaceType, VectorType, DiffusionFactorType, DiffusionTensorType, GridViewType> + LocalNonconformityEstimator; + typedef LocalResidualESV2007<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, GridViewType> + LocalResidualEstimator; + typedef LocalDiffusiveFluxESV2007<SpaceType, VectorType, DiffusionFactorType, DiffusionTensorType, GridViewType> + LocalDiffusiveFluxEstimator; public: - static const bool available = true; - - typedef typename ProblemType::RangeFieldType RangeFieldType; + using typename BaseType::EntityType; + using typename BaseType::ReturnType; - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem) + static std::string id() { - LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, GridType> eta_nc(space, vector, problem); - LocalResidualESV2007<SpaceType, VectorType, ProblemType, GridType> eta_r(space, problem); - LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, GridType> eta_df(space, vector, problem); - eta_nc.prepare(); - eta_r.prepare(); - eta_df.prepare(); + return ESV2007_id(); + } - RangeFieldType eta_squared(0.0); + static ReturnType estimate(const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), space, vector, force, diffusion_factor, diffusion_tensor, over_int); + } - const auto& grid_view = space.grid_view(); - const auto entity_it_end = grid_view.template end<0>(); - for (auto entity_it = grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { - const auto& entity = *entity_it; - eta_squared += - eta_nc.compute_locally(entity) - + std::pow(std::sqrt(eta_r.compute_locally(entity)) + std::sqrt(eta_df.compute_locally(entity)), 2); - } - return std::sqrt(eta_squared); - } // ... estimate(...) + static ReturnType estimate(const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); + } - static Stuff::LA::CommonDenseVector<RangeFieldType> estimate_local(const SpaceType& space, const VectorType& vector, - const ProblemType& problem) + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) { - LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, GridType> eta_nc(space, vector, problem); - LocalResidualESV2007<SpaceType, VectorType, ProblemType, GridType> eta_r(space, problem); - LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, GridType> eta_df(space, vector, problem); - eta_nc.prepare(); - eta_r.prepare(); - eta_df.prepare(); + return estimate(grid_view, space, vector, force, diffusion_factor, diffusion_factor, diffusion_tensor, over_int); + } - const auto& grid_view = space.grid_view(); - Stuff::LA::CommonDenseVector<RangeFieldType> local_indicators( - boost::numeric_cast<size_t>(grid_view.indexSet().size(0)), 0.0); - RangeFieldType eta_squared = 0.0; + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + ThisType estimator(grid_view, + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); + Stuff::Grid::Walker<GridViewType> grid_walker(grid_view); + grid_walker.add(estimator); + grid_walker.walk(); + return std::sqrt(estimator.result()); + } // ... estimate(...) + static Stuff::LA::CommonDenseVector<ReturnType> + estimate_local(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + Stuff::LA::CommonDenseVector<ReturnType> local_indicators(boost::numeric_cast<size_t>(grid_view.indexSet().size(0)), + 0.0); + ReturnType eta_squared = 0.0; + ThisType estimator(grid_view, + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); const auto entity_it_end = grid_view.template end<0>(); for (auto entity_it = grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { - const auto& entity = *entity_it; - const auto index = grid_view.indexSet().index(entity); - const RangeFieldType eta_t_squared = - eta_nc.compute_locally(entity) - + std::pow(std::sqrt(eta_r.compute_locally(entity)) + std::sqrt(eta_df.compute_locally(entity)), 2); - local_indicators[index] = eta_t_squared; + const auto& entity = *entity_it; + const auto index = grid_view.indexSet().index(entity); + const auto eta_t_squared = estimator.compute_locally(entity); + local_indicators[index] = eta_t_squared; eta_squared += eta_t_squared; } for (auto& element : local_indicators) element /= eta_squared; return local_indicators; } // ... estimate_local(...) -}; // class ESV2007< ..., ALUGrid< 2, 2, simplex, conforming >, ... > -#endif // HAVE_ALUGRID + ESV2007(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor_norm, const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + : eta_nc_(grid_view, space, vector, diffusion_factor_norm, diffusion_tensor, over_int) + , eta_r_(grid_view, space, vector, force, diffusion_factor_norm, diffusion_tensor, over_int) + , eta_df_(grid_view, space, vector, diffusion_factor_norm, diffusion_factor_reconstruction, diffusion_tensor, + over_int) + , result_(0.0) + { + } + virtual void prepare() override final + { + eta_nc_.prepare(); + eta_r_.prepare(); + eta_df_.prepare(); + } -class ESV2007AlternativeSummationBase -{ -public: - static std::string id() + virtual ReturnType compute_locally(const EntityType& entity) override { - return "eta_ESV2007_alt"; + return eta_nc_.compute_locally(entity) + + std::pow(std::sqrt(eta_r_.compute_locally(entity)) + std::sqrt(eta_df_.compute_locally(entity)), 2); } -}; + virtual void apply_local(const EntityType& entity) override final + { + result_ += compute_locally(entity); + } -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class ESV2007AlternativeSummation : public ESV2007AlternativeSummationBase -{ -public: - static const bool available = false; -}; + virtual ReturnType result() const override + { + return result_; + } +protected: + LocalNonconformityEstimator eta_nc_; + LocalResidualEstimator eta_r_; + LocalDiffusiveFluxEstimator eta_df_; + ReturnType result_; +}; // class ESV2007 -#if HAVE_ALUGRID -template <class SpaceType, class VectorType, class ProblemType> -class ESV2007AlternativeSummation<SpaceType, VectorType, ProblemType, ALUGrid<2, 2, simplex, conforming>> - : public ESV2007AlternativeSummationBase +template <class SpaceType, class VectorType, class ForceType, class DiffusionFactorType, class DiffusionTensorType, + class GridViewType = typename SpaceType::GridViewType> +class ESV2007AlternativeSummation + : public ESV2007<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, GridViewType> { - typedef ALUGrid<2, 2, simplex, conforming> GridType; + typedef ESV2007AlternativeSummation<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, + GridViewType> ThisType; + typedef ESV2007<SpaceType, VectorType, ForceType, DiffusionFactorType, DiffusionTensorType, GridViewType> BaseType; public: - static const bool available = true; - - typedef typename ProblemType::RangeFieldType RangeFieldType; + using typename BaseType::ReturnType; + using typename BaseType::EntityType; - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem) + static std::string id() { - LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, GridType> eta_nc(space, vector, problem); - LocalResidualESV2007<SpaceType, VectorType, ProblemType, GridType> eta_r(space, problem); - LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, GridType> eta_df(space, vector, problem); - eta_nc.prepare(); - eta_r.prepare(); - eta_df.prepare(); + return ESV2007_alternative_summation_id(); + } - RangeFieldType eta_nc_squared(0.0); - RangeFieldType eta_r_squared(0.0); - RangeFieldType eta_df_squared(0.0); + static ReturnType estimate(const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor, + const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), space, vector, force, diffusion_factor, diffusion_tensor, over_int); + } - const auto& grid_view = space.grid_view(); - const auto entity_it_end = grid_view.template end<0>(); - for (auto entity_it = grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { - const auto& entity = *entity_it; - eta_nc_squared += eta_nc.compute_locally(entity); - eta_r_squared += eta_r.compute_locally(entity); - eta_df_squared += eta_df.compute_locally(entity); - } - return std::sqrt(eta_nc_squared) + std::sqrt(eta_r_squared) + std::sqrt(eta_df_squared); - } // ... estimate(...) + static ReturnType estimate(const SpaceType& space, const VectorType& vector, const ForceType& force, + const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + return estimate(space.grid_view(), + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); + } - static Stuff::LA::CommonDenseVector<RangeFieldType> estimate_local(const SpaceType& space, const VectorType& vector, - const ProblemType& problem) + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) { - LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, GridType> eta_nc(space, vector, problem); - LocalResidualESV2007<SpaceType, VectorType, ProblemType, GridType> eta_r(space, problem); - LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, GridType> eta_df(space, vector, problem); - eta_nc.prepare(); - eta_r.prepare(); - eta_df.prepare(); + return estimate(grid_view, space, vector, force, diffusion_factor, diffusion_factor, diffusion_tensor, over_int); + } - const auto grid_view = space.grid_view(); - Stuff::LA::CommonDenseVector<RangeFieldType> local_indicators( - boost::numeric_cast<size_t>(grid_view.indexSet().size(0)), 0.0); - RangeFieldType eta_nc_squared(0.0); - RangeFieldType eta_r_squared(0.0); - RangeFieldType eta_df_squared(0.0); + static ReturnType estimate(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + ThisType estimator(grid_view, + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); + Stuff::Grid::Walker<GridViewType> grid_walker(grid_view); + grid_walker.add(estimator); + grid_walker.walk(); + return std::sqrt(estimator.result()); + } // ... estimate(...) + static Stuff::LA::CommonDenseVector<ReturnType> + estimate_local(const GridViewType& grid_view, const SpaceType& space, const VectorType& vector, + const ForceType& force, const DiffusionFactorType& diffusion_factor_norm, + const DiffusionFactorType& diffusion_factor_reconstruction, + const DiffusionTensorType& diffusion_tensor, const size_t over_int = over_integrate) + { + Stuff::LA::CommonDenseVector<ReturnType> local_indicators(boost::numeric_cast<size_t>(grid_view.indexSet().size(0)), + 0.0); + ThisType estimator(grid_view, + space, + vector, + force, + diffusion_factor_norm, + diffusion_factor_reconstruction, + diffusion_tensor, + over_int); const auto entity_it_end = grid_view.template end<0>(); for (auto entity_it = grid_view.template begin<0>(); entity_it != entity_it_end; ++entity_it) { - const auto& entity = *entity_it; - const auto index = grid_view.indexSet().index(entity); - const RangeFieldType eta_nc_t_squared = eta_nc.compute_locally(entity); - const RangeFieldType eta_r_t_squared = eta_r.compute_locally(entity); - const RangeFieldType eta_df_t_squared = eta_df.compute_locally(entity); - eta_nc_squared += eta_nc_t_squared; - eta_r_squared += eta_r_t_squared; - eta_df_squared += eta_df_t_squared; - local_indicators[index] = 3.0 * (eta_nc_t_squared + eta_r_t_squared + eta_df_t_squared); + const auto& entity = *entity_it; + const auto index = grid_view.indexSet().index(entity); + local_indicators[index] = 3.0 * estimator.compute_locally(entity); } - const RangeFieldType eta_squared = - std::pow(std::sqrt(eta_nc_squared) + std::sqrt(eta_r_squared) + std::sqrt(eta_df_squared), 2); for (auto& element : local_indicators) - element /= eta_squared; + element /= estimator.result(); return local_indicators; } // ... estimate_local(...) -}; // class ESV2007AlternativeSummation< ..., ALUGrid< 2, 2, simplex, conforming >, ... > - -#endif // HAVE_ALUGRID - - -} // namespace SWIPDG -} // namespace internal - - -template <class SpaceType, class VectorType, class ProblemType, class GridType> -class SWIPDG -{ -public: - typedef typename ProblemType::RangeFieldType RangeFieldType; -private: - template <class IndividualEstimator, bool available = false> - class Caller - { - public: - static std::vector<std::string> append(std::vector<std::string> in) - { - return in; - } - - static bool equals(const std::string& /*type*/) - { - return false; - } - - static RangeFieldType estimate(const SpaceType& /*space*/, const VectorType& /*vector*/, - const ProblemType& /*problem*/) - { - DUNE_THROW(Stuff::Exceptions::internal_error, "This should not happen!"); - return RangeFieldType(0); - } - - static Stuff::LA::CommonDenseVector<RangeFieldType> - estimate_local(const SpaceType& /*space*/, const VectorType& /*vector*/, const ProblemType& /*problem*/) - { - DUNE_THROW(Stuff::Exceptions::internal_error, "This should not happen!"); - return Stuff::LA::CommonDenseVector<RangeFieldType>(); - } - }; // class Caller - - template <class IndividualEstimator> - class Caller<IndividualEstimator, true> + template <class... Args> + ESV2007AlternativeSummation(Args&&... args) + : BaseType(std::forward<Args>(args)...) + , eta_nc_squared_(0.0) + , eta_r_squared_(0.0) + , eta_df_squared_(0.0) { - public: - static std::vector<std::string> append(std::vector<std::string> in) - { - in.push_back(IndividualEstimator::id()); - return in; - } - - static bool equals(const std::string& type) - { - return IndividualEstimator::id() == type; - } - - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem) - { - return IndividualEstimator::estimate(space, vector, problem); - } - - static Stuff::LA::CommonDenseVector<RangeFieldType> estimate_local(const SpaceType& space, const VectorType& vector, - const ProblemType& problem) - { - return IndividualEstimator::estimate_local(space, vector, problem); - } - }; // class Caller< ..., true > - - template <class IndividualEstimator> - static std::vector<std::string> call_append(std::vector<std::string> in) - { - return Caller<IndividualEstimator, IndividualEstimator::available>::append(in); - } - - template <class IndividualEstimator> - static bool call_equals(const std::string& type) - { - return Caller<IndividualEstimator, IndividualEstimator::available>::equals(type); } - template <class IndividualEstimator> - static RangeFieldType call_estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem) + virtual ReturnType compute_locally(const EntityType& entity) override final { - return Caller<IndividualEstimator, IndividualEstimator::available>::estimate(space, vector, problem); + const auto eta_nc_t_squared = eta_nc_.compute_locally(entity); + const auto eta_r_t_squared = eta_r_.compute_locally(entity); + const auto eta_df_t_squared = eta_df_.compute_locally(entity); + eta_nc_squared_ += eta_nc_t_squared; + eta_r_squared_ += eta_r_t_squared; + eta_df_squared_ += eta_df_t_squared; + return eta_nc_t_squared + eta_r_t_squared + eta_df_t_squared; } - template <class IndividualEstimator> - static Stuff::LA::CommonDenseVector<RangeFieldType> - call_estimate_local(const SpaceType& space, const VectorType& vector, const ProblemType& problem) + virtual ReturnType result() const override { - return Caller<IndividualEstimator, IndividualEstimator::available>::estimate_local(space, vector, problem); + return std::sqrt(eta_nc_squared_) + std::sqrt(eta_r_squared_) + std::sqrt(eta_df_squared_); } - typedef internal::SWIPDG::LocalNonconformityESV2007<SpaceType, VectorType, ProblemType, GridType> - LocalNonconformityESV2007Type; - typedef internal::SWIPDG::LocalResidualESV2007<SpaceType, VectorType, ProblemType, GridType> LocalResidualESV2007Type; - typedef internal::SWIPDG::LocalResidualESV2007Star<SpaceType, VectorType, ProblemType, GridType> - LocalResidualESV2007StarType; - typedef internal::SWIPDG::LocalDiffusiveFluxESV2007<SpaceType, VectorType, ProblemType, GridType> - LocalDiffusiveFluxESV2007Type; - typedef internal::SWIPDG::ESV2007<SpaceType, VectorType, ProblemType, GridType> ESV2007Type; - typedef internal::SWIPDG::ESV2007AlternativeSummation<SpaceType, VectorType, ProblemType, GridType> - ESV2007AlternativeSummationType; - -public: - static std::vector<std::string> available() - { - std::vector<std::string> tmp; - tmp = call_append<LocalNonconformityESV2007Type>(tmp); - tmp = call_append<LocalResidualESV2007Type>(tmp); - tmp = call_append<LocalResidualESV2007StarType>(tmp); - tmp = call_append<LocalDiffusiveFluxESV2007Type>(tmp); - tmp = call_append<ESV2007Type>(tmp); - tmp = call_append<ESV2007AlternativeSummationType>(tmp); - return tmp; - } // ... available(...) - - static std::vector<std::string> available_local() - { - std::vector<std::string> tmp; - tmp = call_append<ESV2007Type>(tmp); - tmp = call_append<ESV2007AlternativeSummationType>(tmp); - return tmp; - } // ... available_local(...) - - static RangeFieldType estimate(const SpaceType& space, const VectorType& vector, const ProblemType& problem, - const std::string type) - { - if (call_equals<LocalNonconformityESV2007Type>(type)) - return call_estimate<LocalNonconformityESV2007Type>(space, vector, problem); - else if (call_equals<LocalResidualESV2007Type>(type)) - return call_estimate<LocalResidualESV2007Type>(space, vector, problem); - else if (call_equals<LocalResidualESV2007StarType>(type)) - return call_estimate<LocalResidualESV2007StarType>(space, vector, problem); - else if (call_equals<LocalDiffusiveFluxESV2007Type>(type)) - return call_estimate<LocalDiffusiveFluxESV2007Type>(space, vector, problem); - else if (call_equals<ESV2007Type>(type)) - return call_estimate<ESV2007Type>(space, vector, problem); - else if (call_equals<ESV2007AlternativeSummationType>(type)) - return call_estimate<ESV2007AlternativeSummationType>(space, vector, problem); - else - DUNE_THROW(Stuff::Exceptions::you_are_using_this_wrong, - "Requested type '" << type << "' is not one of available()!"); - } // ... estimate(...) - - static Stuff::LA::CommonDenseVector<RangeFieldType> estimate_local(const SpaceType& space, const VectorType& vector, - const ProblemType& problem, const std::string type) - { - if (call_equals<ESV2007Type>(type)) - return call_estimate_local<ESV2007Type>(space, vector, problem); - else if (call_equals<ESV2007AlternativeSummationType>(type)) - return call_estimate_local<ESV2007AlternativeSummationType>(space, vector, problem); - else - DUNE_THROW(Stuff::Exceptions::you_are_using_this_wrong, - "Requested type '" << type << "' is not one of available_local()!"); - } // ... estimate_local(...) -}; // class SWIPDG - - -} // namespace Discretizations -} // namespace Estimators -} // namespace HDD +private: + using BaseType::eta_nc_; + using BaseType::eta_r_; + using BaseType::eta_df_; + ReturnType eta_nc_squared_; + ReturnType eta_r_squared_; + ReturnType eta_df_squared_; +}; // class ESV2007AlternativeSummation + + +} // namespace SwipdgFluxreconstrutionEstimators +} // namespace LinearElliptic +} // namespace GDT } // namespace Dune -#endif // DUNE_HDD_LINEARELLIPTIC_ESTIMATORS_SWIPDG_HH +#endif // DUNE_GDT_TEST_LINEARELLIPTIC_ESTIMATORS_SWIPDG_FLUXRECONSTRUCTION_HH -- GitLab