diff --git a/dune/gdt/operators/advection-dg.hh b/dune/gdt/operators/advection-dg.hh index 60eb3661e1fb9405f8b8a279fdd7c5f2189e8971..3ec5e91fb6012ee40644422d9a6abfd375fcab29 100644 --- a/dune/gdt/operators/advection-dg.hh +++ b/dune/gdt/operators/advection-dg.hh @@ -39,72 +39,39 @@ namespace GDT { * * \sa OperatorInterface */ -template <class AssemblyGridView, - class SV, - size_t m = 1, - class SF = double, - class SGV = AssemblyGridView, - class RF = SF, - class RGV = AssemblyGridView, - class RV = SV> -class AdvectionDgOperator : public OperatorInterface<SV, - SGV, - m, - 1, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - m, - 1, - RF, - RGV, - RV> +template <class M, class SGV, size_t m = 1, class RGV = SGV> +class AdvectionDgOperator : public OperatorInterface<M, SGV, m, 1, m, 1, RGV> { - // No need to check the rest, is done in OperatorInterface. - static_assert(XT::Grid::is_view<AssemblyGridView>::value, ""); - static_assert((AssemblyGridView::dimension == SGV::dimension) && (AssemblyGridView::dimension == RGV::dimension), ""); - - using ThisType = AdvectionDgOperator<AssemblyGridView, SV, m, SF, SGV, RF, RGV, RV>; - using BaseType = OperatorInterface<SV, - SGV, - m, - 1, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - m, - 1, - RF, - RGV, - RV>; + using ThisType = AdvectionDgOperator<M, SGV, m, RGV>; + using BaseType = OperatorInterface<M, SGV, m, 1, m, 1, RGV>; protected: - using AGV = AssemblyGridView; - static const constexpr size_t d = AssemblyGridView::dimension; - using D = typename AssemblyGridView::ctype; + static const constexpr size_t d = SGV::dimension; + using D = typename SGV::ctype; public: - using AssemblyGridViewType = AssemblyGridView; - using NumericalFluxType = NumericalFluxInterface<d, m, RF>; + using typename BaseType::F; + using typename BaseType::V; - using I = XT::Grid::extract_intersection_t<AGV>; + using NumericalFluxType = NumericalFluxInterface<d, m, F>; + + using I = XT::Grid::extract_intersection_t<SGV>; using BoundaryTreatmentByCustomNumericalFluxOperatorType = - LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator<I, SV, SGV, m, SF, RF, RGV, RV>; + LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator<I, V, SGV, m, F, F, RGV, V>; using BoundaryTreatmentByCustomExtrapolationOperatorType = - LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator<I, SV, SGV, m, SF, RF, RGV, RV>; + LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator<I, V, SGV, m, F, F, RGV, V>; using typename BaseType::SourceSpaceType; - using typename BaseType::SourceVectorType; - + using typename BaseType::VectorType; using typename BaseType::RangeSpaceType; - using typename BaseType::RangeVectorType; using typename BaseType::RangeFunctionType; AdvectionDgOperator( - const AssemblyGridViewType& assembly_grid_view, + const SGV& assembly_grid_view, const NumericalFluxType& numerical_flux, const SourceSpaceType& source_space, const RangeSpaceType& range_space, - const XT::Grid::IntersectionFilter<AGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<AGV>(), - const size_t jump_indicator_component = 0) + const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>()) : BaseType(numerical_flux.parameter_type()) , assembly_grid_view_(assembly_grid_view) , numerical_flux_(numerical_flux.copy()) @@ -138,7 +105,7 @@ public: append(typename BoundaryTreatmentByCustomNumericalFluxOperatorType::LambdaType numerical_boundary_treatment_flux, const int numerical_boundary_treatment_flux_order, const XT::Common::ParameterType& boundary_treatment_parameter_type = {}, - const XT::Grid::IntersectionFilter<AGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<AGV>()) + const XT::Grid::IntersectionFilter<SGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<SGV>()) { boundary_treatments_by_custom_numerical_flux_.emplace_back( new BoundaryTreatmentByCustomNumericalFluxOperatorType(numerical_boundary_treatment_flux, @@ -150,7 +117,7 @@ public: ThisType& append(typename BoundaryTreatmentByCustomExtrapolationOperatorType::LambdaType extrapolation, const XT::Common::ParameterType& extrapolation_parameter_type = {}, - const XT::Grid::IntersectionFilter<AGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<AGV>()) + const XT::Grid::IntersectionFilter<SGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<SGV>()) { boundary_treatments_by_custom_extrapolation_.emplace_back( new BoundaryTreatmentByCustomExtrapolationOperatorType( @@ -162,19 +129,19 @@ public: /// \} protected: - void append_standard_contributions(LocalizableOperatorBase<AGV, SV, m, 1, SF, SGV, m, 1, RF, RGV, RV>& localizable_op, + void append_standard_contributions(LocalizableOperatorBase<SGV, V, m, 1, F, SGV, m, 1, F, RGV, V>& localizable_op, const XT::Common::Parameter& param) const { // element contributions - localizable_op.append(LocalAdvectionDgVolumeOperator<SV, SGV, m, SF, RF, RGV, RV>(numerical_flux_->flux()), param); + localizable_op.append(LocalAdvectionDgVolumeOperator<V, SGV, m, F, F, RGV, V>(numerical_flux_->flux()), param); // contributions from inner intersections - localizable_op.append(LocalAdvectionDgCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(*numerical_flux_), + localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_), param, - XT::Grid::ApplyOn::InnerIntersectionsOnce<AGV>()); + XT::Grid::ApplyOn::InnerIntersectionsOnce<SGV>()); // contributions from periodic boundaries - localizable_op.append(LocalAdvectionDgCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(*numerical_flux_), + localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_), param, - *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<AGV>() && !(*periodicity_exception_))); + *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<SGV>() && !(*periodicity_exception_))); // contributions from other boundaries by custom numerical flux for (const auto& boundary_treatment : boundary_treatments_by_custom_numerical_flux_) { const auto& boundary_op = *boundary_treatment.first; @@ -189,28 +156,27 @@ protected: } } // ... append_standard_contributions(...) - void append_local_mass_matrix_inversion( - LocalizableOperatorBase<AGV, SV, m, 1, SF, SGV, m, 1, RF, RGV, RV>& localizable_op, - RangeFunctionType& range_function) const + void + append_local_mass_matrix_inversion(LocalizableOperatorBase<SGV, V, m, 1, F, SGV, m, 1, F, RGV, V>& localizable_op, + RangeFunctionType& range_function) const { localizable_op.append([&](const auto& element) { // (creating these objects before the grid walk and reusing them would be more efficient, but not thread safe) auto local_range = range_function.local_discrete_function(element); const auto& range_basis = local_range->basis(); using E = XT::Grid::extract_entity_t<RGV>; - const LocalElementIntegralBilinearForm<E, m, 1, RF, RF> local_l2_bilinear_form( - LocalElementProductIntegrand<E, m, 1, RF, RF>(1.)); + const LocalElementIntegralBilinearForm<E, m, 1, F, F> local_l2_bilinear_form( + LocalElementProductIntegrand<E, m, 1, F, F>(1.)); local_range->dofs().assign_from(XT::LA::solve( - XT::LA::convert_to<XT::LA::CommonDenseMatrix<RF>>(local_l2_bilinear_form.apply2(range_basis, range_basis)), - XT::LA::convert_to<XT::LA::CommonDenseVector<RF>>(local_range->dofs()), - {{"type", XT::LA::SolverOptions<XT::LA::CommonDenseMatrix<RF>>::types().at(0)}, + XT::LA::convert_to<XT::LA::CommonDenseMatrix<F>>(local_l2_bilinear_form.apply2(range_basis, range_basis)), + XT::LA::convert_to<XT::LA::CommonDenseVector<F>>(local_range->dofs()), + {{"type", XT::LA::SolverOptions<XT::LA::CommonDenseMatrix<F>>::types().at(0)}, {"post_check_solves_system", "-1"}})); }); } // ... append_local_mass_matrix_inversion(...) public: - void - apply(const SourceVectorType& source, RangeVectorType& range, const XT::Common::Parameter& param = {}) const override + void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override { // some checks DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!"); @@ -234,31 +200,30 @@ public: } // ... apply(...) protected: - const AssemblyGridViewType& assembly_grid_view_; + const SGV assembly_grid_view_; const std::unique_ptr<const NumericalFluxType> numerical_flux_; const SourceSpaceType& source_space_; const RangeSpaceType& range_space_; - std::unique_ptr<XT::Grid::IntersectionFilter<AssemblyGridViewType>> periodicity_exception_; + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>> periodicity_exception_; std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>, - std::unique_ptr<XT::Grid::IntersectionFilter<AGV>>>> + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>> boundary_treatments_by_custom_numerical_flux_; std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomExtrapolationOperatorType>, - std::unique_ptr<XT::Grid::IntersectionFilter<AGV>>>> + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>> boundary_treatments_by_custom_extrapolation_; }; // class AdvectionDgOperator -template <class VectorType, class AGV, size_t m, class RF, class SGV, class SF, class RGV> -std::enable_if_t<XT::LA::is_vector<VectorType>::value, - AdvectionDgOperator<AGV, VectorType, m, SF, SGV, RF, RGV, VectorType>> +template <class MatrixType, class SGV, size_t m, class F, class RGV> +std::enable_if_t<XT::LA::is_matrix<MatrixType>::value, AdvectionDgOperator<MatrixType, SGV, m, RGV>> make_advection_dg_operator( - const AGV& assembly_grid_view, - const NumericalFluxInterface<AGV::dimension, m, RF>& numerical_flux, - const SpaceInterface<SGV, m, 1, SF>& source_space, - const SpaceInterface<RGV, m, 1, RF>& range_space, - const XT::Grid::IntersectionFilter<AGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<AGV>()) + const SGV& assembly_grid_view, + const NumericalFluxInterface<SGV::dimension, m, F>& numerical_flux, + const SpaceInterface<SGV, m, 1, F>& source_space, + const SpaceInterface<RGV, m, 1, F>& range_space, + const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>()) { - return AdvectionDgOperator<AGV, VectorType, m, SF, SGV, RF, RGV, VectorType>( + return AdvectionDgOperator<MatrixType, SGV, m, RGV>( assembly_grid_view, numerical_flux, source_space, range_space, periodicity_exception); } @@ -269,37 +234,27 @@ make_advection_dg_operator( * \sa AdvectionDgOperator * \sa [FK2007, Sec. 6.2] for nu_1 */ -template <class AGV, - class SV, - size_t m = 1, - class SF = double, - class SGV = AGV, - class RF = SF, - class RGV = AGV, - class RV = SV> -class AdvectionDgArtificialViscosityOperator : public AdvectionDgOperator<AGV, SV, m, SF, SGV, RF, RGV, RV> +template <class M, class SGV, size_t m = 1, class RGV = SGV> +class AdvectionDgArtificialViscosityOperator : public AdvectionDgOperator<M, SGV, m, RGV> { - using ThisType = AdvectionDgArtificialViscosityOperator<AGV, SV, m, SF, SGV, RF, RGV, RV>; - using BaseType = AdvectionDgOperator<AGV, SV, m, SF, SGV, RF, RGV, RV>; + using ThisType = AdvectionDgArtificialViscosityOperator<M, SGV, m, RGV>; + using BaseType = AdvectionDgOperator<M, SGV, m, RGV>; using BaseType::d; using typename BaseType::D; public: - using typename BaseType::AssemblyGridViewType; using typename BaseType::NumericalFluxType; - using typename BaseType::SourceSpaceType; using typename BaseType::RangeSpaceType; - - using typename BaseType::SourceVectorType; - using typename BaseType::RangeVectorType; + using typename BaseType::VectorType; + using typename BaseType::F; AdvectionDgArtificialViscosityOperator( - const AssemblyGridViewType& assembly_grid_view, + const SGV& assembly_grid_view, const NumericalFluxType& numerical_flux, const SourceSpaceType& source_space, const RangeSpaceType& range_space, - const XT::Grid::IntersectionFilter<AGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<AGV>(), + const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>(), const double nu_1 = 0.2, const double alpha_1 = 1., const size_t jump_indicator_component = 0) @@ -314,9 +269,7 @@ public: using BaseType::apply; - void apply(const SourceVectorType& source, - RangeVectorType& range, - const XT::Common::Parameter& param = {}) const override final + void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override final { // some checks DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!"); @@ -331,7 +284,7 @@ public: auto walker = XT::Grid::make_walker(assembly_grid_view_); walker.append([&](const auto& element) { auto local_source = source_function.local_discrete_function(element); - using E = XT::Grid::extract_entity_t<AGV>; + using E = XT::Grid::extract_entity_t<SGV>; const auto average = LocalElementIntegralFunctional<E>(LocalElementIdentityIntegrand<E>()).apply(*local_source)[0] / element.geometry().volume(); local_source->dofs().assign_from(source_space_.finite_element(element.geometry().type()) @@ -348,7 +301,7 @@ public: // compute jump indicators for artificial viscosity at detected shocks, see DF2015, p. 449, (8.180) // we need a thread safe vector with one entry per grid element and just use a scalar FV discrete function const auto fv_space = make_finite_volume_space<1>(this->assembly_grid_view_); - auto jump_indicators = make_discrete_function<XT::LA::CommonDenseVector<RF>>(fv_space); + auto jump_indicators = make_discrete_function<XT::LA::CommonDenseVector<F>>(fv_space); localizable_op.append([&](const auto& element) { // compute jump indicator (8.176) double element_jump_indicator = 0; diff --git a/dune/gdt/operators/advection-fv.hh b/dune/gdt/operators/advection-fv.hh index 0fd409d5ea8e05eb86abef5ea23891c3c21d748e..e6e8655c78cadacade9698adcc853a50c61e9a2d 100644 --- a/dune/gdt/operators/advection-fv.hh +++ b/dune/gdt/operators/advection-fv.hh @@ -28,66 +28,34 @@ namespace GDT { * * \sa OperatorInterface */ -template <class AssemblyGridView, - class SV, - size_t m = 1, - class SF = double, - class SGV = AssemblyGridView, - class RF = SF, - class RGV = AssemblyGridView, - class RV = SV> -class AdvectionFvOperator : public OperatorInterface<SV, - SGV, - m, - 1, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - m, - 1, - RF, - RGV, - RV> +template <class M, class SGV, size_t m = 1, class RGV = SGV> +class AdvectionFvOperator : public OperatorInterface<M, SGV, m, 1, m, 1, RGV> { - // No need to check the rest, is done in OperatorInterface. - static_assert(XT::Grid::is_view<AssemblyGridView>::value, ""); - static_assert((AssemblyGridView::dimension == SGV::dimension) && (AssemblyGridView::dimension == RGV::dimension), ""); - - using ThisType = AdvectionFvOperator<AssemblyGridView, SV, m, SF, SGV, RF, RGV, RV>; - using BaseType = OperatorInterface<SV, - SGV, - m, - 1, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - m, - 1, - RF, - RGV, - RV>; - using AGV = AssemblyGridView; + using ThisType = AdvectionFvOperator<M, SGV, m, RGV>; + using BaseType = OperatorInterface<M, SGV, m, 1, m, 1, RGV>; public: - using AssemblyGridViewType = AssemblyGridView; - using NumericalFluxType = NumericalFluxInterface<AssemblyGridView::dimension, m, RF>; + using typename BaseType::V; + using typename BaseType::F; - using I = XT::Grid::extract_intersection_t<AGV>; + using NumericalFluxType = NumericalFluxInterface<SGV::dimension, m, F>; + + using I = XT::Grid::extract_intersection_t<SGV>; using BoundaryTreatmentByCustomNumericalFluxOperatorType = - LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator<I, SV, SGV, m, SF, RF, RGV, RV>; + LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator<I, V, SGV, m, F, F, RGV, V>; using BoundaryTreatmentByCustomExtrapolationOperatorType = - LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, SV, SGV, m, SF, RF, RGV, RV>; + LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, V, SGV, m, F, F, RGV, V>; using typename BaseType::SourceSpaceType; using typename BaseType::RangeSpaceType; - - using typename BaseType::SourceVectorType; - using typename BaseType::RangeVectorType; + using typename BaseType::VectorType; AdvectionFvOperator( - const AssemblyGridViewType& assembly_grid_view, + const SGV& assembly_grid_view, const NumericalFluxType& numerical_flux, const SourceSpaceType& source_space, const RangeSpaceType& range_space, - const XT::Grid::IntersectionFilter<AGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<AGV>()) + const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>()) : BaseType(numerical_flux.parameter_type()) , assembly_grid_view_(assembly_grid_view) , numerical_flux_(numerical_flux.copy()) @@ -120,7 +88,7 @@ public: ThisType& append(typename BoundaryTreatmentByCustomNumericalFluxOperatorType::LambdaType numerical_boundary_treatment_flux, const XT::Common::ParameterType& boundary_treatment_parameter_type = {}, - const XT::Grid::IntersectionFilter<AGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<AGV>()) + const XT::Grid::IntersectionFilter<SGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<SGV>()) { boundary_treatments_by_custom_numerical_flux_.emplace_back( new BoundaryTreatmentByCustomNumericalFluxOperatorType(numerical_boundary_treatment_flux, @@ -131,7 +99,7 @@ public: ThisType& append(typename BoundaryTreatmentByCustomExtrapolationOperatorType::LambdaType extrapolation, const XT::Common::ParameterType& extrapolation_parameter_type = {}, - const XT::Grid::IntersectionFilter<AGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<AGV>()) + const XT::Grid::IntersectionFilter<SGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<SGV>()) { boundary_treatments_by_custom_extrapolation_.emplace_back( new BoundaryTreatmentByCustomExtrapolationOperatorType( @@ -144,9 +112,7 @@ public: using BaseType::apply; - void apply(const SourceVectorType& source, - RangeVectorType& range, - const XT::Common::Parameter& param = {}) const override final + void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override final { // some checks DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!"); @@ -159,13 +125,13 @@ public: // set up the actual operator auto localizable_op = make_localizable_operator(assembly_grid_view_, source_function, range_function); // contributions from inner intersections - localizable_op.append(LocalAdvectionFvCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(*numerical_flux_), + localizable_op.append(LocalAdvectionFvCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_), param, - XT::Grid::ApplyOn::InnerIntersectionsOnce<AGV>()); + XT::Grid::ApplyOn::InnerIntersectionsOnce<SGV>()); // contributions from periodic boundaries - localizable_op.append(LocalAdvectionFvCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(*numerical_flux_), + localizable_op.append(LocalAdvectionFvCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_), param, - *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<AGV>() && !(*periodicity_exception_))); + *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<SGV>() && !(*periodicity_exception_))); // contributions from other boundaries by custom numerical flux for (const auto& boundary_treatment : boundary_treatments_by_custom_numerical_flux_) { const auto& boundary_op = *boundary_treatment.first; @@ -184,31 +150,30 @@ public: } // ... apply(...) private: - const AssemblyGridViewType& assembly_grid_view_; + const SGV assembly_grid_view_; const std::unique_ptr<const NumericalFluxType> numerical_flux_; const SourceSpaceType& source_space_; const RangeSpaceType& range_space_; - std::unique_ptr<XT::Grid::IntersectionFilter<AssemblyGridViewType>> periodicity_exception_; + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>> periodicity_exception_; std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>, - std::unique_ptr<XT::Grid::IntersectionFilter<AGV>>>> + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>> boundary_treatments_by_custom_numerical_flux_; std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomExtrapolationOperatorType>, - std::unique_ptr<XT::Grid::IntersectionFilter<AGV>>>> + std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>> boundary_treatments_by_custom_extrapolation_; }; // class AdvectionFvOperator -template <class VectorType, class AGV, size_t m, class RF, class SGV, class SF, class RGV> -std::enable_if_t<XT::LA::is_vector<VectorType>::value, - AdvectionFvOperator<AGV, VectorType, m, SF, SGV, RF, RGV, VectorType>> +template <class MatrixType, class SGV, size_t m, class F, class RGV> +std::enable_if_t<XT::LA::is_matrix<MatrixType>::value, AdvectionFvOperator<MatrixType, SGV, m, RGV>> make_advection_fv_operator( - const AGV& assembly_grid_view, - const NumericalFluxInterface<AGV::dimension, m, RF>& numerical_flux, - const SpaceInterface<SGV, m, 1, SF>& source_space, - const SpaceInterface<RGV, m, 1, RF>& range_space, - const XT::Grid::IntersectionFilter<AGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<AGV>()) + const SGV& assembly_grid_view, + const NumericalFluxInterface<SGV::dimension, m, F>& numerical_flux, + const SpaceInterface<SGV, m, 1, F>& source_space, + const SpaceInterface<RGV, m, 1, F>& range_space, + const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>()) { - return AdvectionFvOperator<AGV, VectorType, m, SF, SGV, RF, RGV, VectorType>( + return AdvectionFvOperator<MatrixType, SGV, m, RGV>( assembly_grid_view, numerical_flux, source_space, range_space, periodicity_exception); } diff --git a/dune/gdt/operators/interfaces.hh b/dune/gdt/operators/interfaces.hh index 568d1aa39364f0101a72da6abc2fa06178d02bc8..201d73b94e1ab6378d59d29e20e0a10447ae102f 100644 --- a/dune/gdt/operators/interfaces.hh +++ b/dune/gdt/operators/interfaces.hh @@ -1,6 +1,6 @@ // This file is part of the dune-gdt project: // https://github.com/dune-community/dune-gdt -// Copyright 2010-2018 dune-gdt developers and contributors. All rights reserved. +// Copyright 2010-2018 dune-gdt developers and contributors. All rights reseVed. // License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) // or GPL-2.0+ (http://opensource.org/licenses/gpl-license) // with "runtime exception" (http://www.dune-project.org/license.html) @@ -19,90 +19,122 @@ #include <dune/xt/common/configuration.hh> #include <dune/xt/common/type_traits.hh> #include <dune/xt/la/container/vector-interface.hh> +#include <dune/xt/la/type_traits.hh> #include <dune/gdt/discretefunction/default.hh> #include <dune/gdt/exceptions.hh> #include <dune/gdt/spaces/interface.hh> +#include <dune/gdt/tools/sparsity-pattern.hh> namespace Dune { namespace GDT { +namespace internal { + + +template <class Matrix, + class SourceGridView, + size_t source_dim, + size_t source_dim_cols, + size_t range_dim, + size_t range_dim_cols, + class RangeGridView> +class AssertArgumentsOfOperatorInterface +{ + static_assert(XT::LA::is_matrix<Matrix>::value, ""); + static_assert(XT::Grid::is_view<SourceGridView>::value, ""); + static_assert(source_dim > 0, ""); + static_assert(source_dim_cols > 0, ""); + static_assert(range_dim > 0, ""); + static_assert(range_dim_cols > 0, ""); + static_assert(XT::Grid::is_view<RangeGridView>::value, ""); + static_assert(SourceGridView::dimension == RangeGridView::dimension, ""); +}; // class AssertArgumentsOfOperatorInterface + + +} // namespace internal + + +// forward, required for the jacobian +template <class M, class SGV, size_t s_r, size_t s_rC, size_t r_r, size_t r_rC, class RGV> +class MatrixOperator; /** * \brief Interface for operators (and two-forms). * - * Considering (discrete) spaces V_h and W_h, and a field K, this interface models + * Considering (discrete) spaces V_h and W_h, and a field F, this interface models * * - operators A: V_h -> W_h and * - * - two-forms B: W_h x V_h -> K. + * - two-forms B: W_h x V_h -> F. * * The source of the operator V_h is the (discrete) space * - * V_h := \{ v_h: \tau_h^s -> SF^{s_r \times s_rC} \} + * V_h := \{ v_h: \tau_h^s -> F^{s_r \times s_rC} \} * * (modelled by SpaceInterface) of functions mapping from a partition \tau_h^s (modelled by SourceGridView) of a - * physical domain to a (possibly vector- or matrix-valued) vector space SF^{s_r \times s_rC} (modelled by SourceField, + * physical domain to a (possibly vector- or matrix-valued) vector space F^{s_r \times s_rC} (modelled by Field, * source_dim and source_dim_cols). The range of the operator W_h (identified with its dual, since we are in the * discrete setting) is the (discrete) space * - * W_h := \{ w_h: \tau_h^r -> RF^{r_r \times r_rC} \} + * W_h := \{ w_h: \tau_h^r -> F^{r_r \times r_rC} \} * * (modelled by SpaceInterface) of functions mapping from a partition \tau_h^r (modelled by RangeGridView) of a - * physical domain to a (possibly vector- or matrix-valued) vector space RF^{r_r \times r_rC} (modelled by RangeField, + * physical domain to a (possibly vector- or matrix-valued) vector space F^{r_r \times r_rC} (modelled by Field, * range_dim and range_dim_cols). * * The functions v_h \in V_h (and w_h \in W_h), to which the operator can be applied to, are represented by their DoF - * vectors (an appropriate vector type derived from XT::LA::VectorInterface modelled by SourceVector (RangeVector)), - * which are then interpreted as discrete functions in the respective source_space or range_space. + * vectors, which are then interpreted as discrete functions in the respective source_space or range_space. * - * The field K of the interpretation of the operator as a two-form (see for instance the default implementation of - * apply2()) is modelled by Field. + * The appropriate vector type (derived from XT::LA::VectorInterface) is automatically deduced from the given + * matrix type (derived from XT::LA::MatrixInterface, modelled by Matrix), as well as the underlying field. + * + * \note In general, one would like to have differente fields for the source vector, the range vector, the matrix and + * the result of apply2(). However, this is postponed in favor of fewer template arguments, until we require it. */ -template <class SourceVector, +template <class Matrix, class SourceGridView, size_t source_dim = 1, size_t source_dim_cols = 1, - class SourceField = double, - class Field = double, size_t range_dim = source_dim, size_t range_dim_cols = source_dim_cols, - class RangeField = double, - class RangeGridView = SourceGridView, - class RangeVector = SourceVector> -class OperatorInterface : public XT::Common::ParametricInterface + class RangeGridView = SourceGridView> +class OperatorInterface : internal::AssertArgumentsOfOperatorInterface<Matrix, + SourceGridView, + source_dim, + source_dim_cols, + range_dim, + range_dim_cols, + RangeGridView>, + public XT::Common::ParametricInterface { - static_assert(XT::LA::is_vector<SourceVector>::value, ""); - static_assert(XT::LA::is_vector<RangeVector>::value, ""); - public: - using SV = SourceVector; + using MatrixType = Matrix; + using VectorType = XT::LA::vector_t<MatrixType>; + using FieldType = typename MatrixType::ScalarType; + + using M = Matrix; + using V = VectorType; + using F = FieldType; + using SGV = SourceGridView; static const constexpr size_t s_r = source_dim; static const constexpr size_t s_rC = source_dim_cols; - using SF = SourceField; - using RV = RangeVector; using RGV = RangeGridView; static const constexpr size_t r_r = range_dim; static const constexpr size_t r_rC = range_dim_cols; - using RF = RangeField; - - using F = Field; - using SourceSpaceType = SpaceInterface<SGV, s_r, s_rC, SF>; - using SourceVectorType = SourceVector; - using SourceFunctionType = DiscreteFunction<SourceVector, SGV, s_r, s_rC, SF>; - using ConstSourceFunctionType = ConstDiscreteFunction<SourceVector, SGV, s_r, s_rC, SF>; + using SourceSpaceType = SpaceInterface<SGV, s_r, s_rC, F>; + using SourceFunctionType = DiscreteFunction<V, SGV, s_r, s_rC, F>; + using ConstSourceFunctionType = ConstDiscreteFunction<V, SGV, s_r, s_rC, F>; - using RangeSpaceType = SpaceInterface<RGV, r_r, r_rC, RF>; - using RangeVectorType = RangeVector; - using RangeFunctionType = DiscreteFunction<RangeVector, RGV, r_r, r_rC, RF>; - using ConstRangeFunctionType = ConstDiscreteFunction<RangeVector, RGV, r_r, r_rC, RF>; + using RangeSpaceType = SpaceInterface<RGV, r_r, r_rC, F>; + using RangeFunctionType = DiscreteFunction<V, RGV, r_r, r_rC, F>; + using ConstRangeFunctionType = ConstDiscreteFunction<V, RGV, r_r, r_rC, F>; - using FieldType = Field; - - using ThisType = OperatorInterface<SV, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, RV>; + using ThisType = OperatorInterface<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; + using MatrixOperatorType = MatrixOperator<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; explicit OperatorInterface(const XT::Common::ParameterType& param_type = {}) : XT::Common::ParametricInterface(param_type) @@ -134,9 +166,8 @@ public: /// \name These methods should be implemented and define the functionality of the operator. /// \{ - virtual void apply(const SourceVectorType& /*source*/, - RangeVectorType& /*range*/, - const XT::Common::Parameter& /*param*/ = {}) const + virtual void + apply(const VectorType& /*source*/, VectorType& /*range*/, const XT::Common::Parameter& /*param*/ = {}) const { DUNE_THROW(Exceptions::operator_error, "This operator cannot be applied!"); } @@ -169,8 +200,8 @@ invert_options(some_type).get<std::string>("type") == some_type return XT::Common::Configuration(); } - virtual void apply_inverse(const RangeVectorType& /*range*/, - SourceVectorType& /*source*/, + virtual void apply_inverse(const VectorType& /*range*/, + VectorType& /*source*/, const XT::Common::Configuration& /*opts*/, const XT::Common::Parameter& /*param*/ = {}) const { @@ -203,21 +234,27 @@ invert_options(some_type).get<std::string>("type") == some_type return XT::Common::Configuration(); } - virtual std::shared_ptr<ThisType> jacobian(const SourceVectorType& /*source*/, - const XT::Common::Configuration& /*opts*/, - const XT::Common::Parameter& /*param*/ = {}) const + /** + * Either appends suitable functors to the jacobian_op (such that the jacobian of this operator is assembled + * additively into jacobian_op) or adds the jacobian of this operator to jacobian_op.matrix(). + * + * \note That you need to call jacobian_op.assemble() to be sure to have the jacobian fully assembled. + **/ + virtual void jacobian(const VectorType& /*source*/, + MatrixOperatorType& /*jacobian_op*/, + const XT::Common::Configuration& /*opts*/, + const XT::Common::Parameter& /*param*/ = {}) const { DUNE_THROW(Exceptions::operator_error, "This operator does not provide a jacobian!"); - return nullptr; } /// \} /// \name These apply variants are provided for convenience. /// \{ - virtual RangeVectorType apply(const SourceVectorType& source, const XT::Common::Parameter& param = {}) const + virtual VectorType apply(const VectorType& source, const XT::Common::Parameter& param = {}) const { - RangeVectorType range(this->range_space().mapper().size(), 0); + VectorType range(this->range_space().mapper().size(), 0); this->apply(source, range, param); return range; } @@ -227,7 +264,7 @@ invert_options(some_type).get<std::string>("type") == some_type /// \{ virtual FieldType - apply2(const SourceVectorType& range, const RangeVectorType& source, const XT::Common::Parameter& param = {}) const + apply2(const VectorType& range, const VectorType& source, const XT::Common::Parameter& param = {}) const { return range.dot(this->apply(source, param)); } @@ -236,40 +273,39 @@ invert_options(some_type).get<std::string>("type") == some_type /// \name These apply_invers variants are provided for convenience. /// \{ - virtual void apply_inverse(const RangeVectorType& range, - SourceVectorType& source, + virtual void apply_inverse(const VectorType& range, + VectorType& source, const std::string& type, const XT::Common::Parameter& param = {}) const { this->apply_inverse(range, source, this->invert_options(type), param); } - virtual void - apply_inverse(const RangeVectorType& range, SourceVectorType& source, const XT::Common::Parameter& param = {}) const + virtual void apply_inverse(const VectorType& range, VectorType& source, const XT::Common::Parameter& param = {}) const { this->apply_inverse(range, source, this->invert_options().at(0), param); } - virtual SourceVectorType apply_inverse(const RangeVectorType& range, - const XT::Common::Configuration& opts, - const XT::Common::Parameter& param = {}) const + virtual VectorType apply_inverse(const VectorType& range, + const XT::Common::Configuration& opts, + const XT::Common::Parameter& param = {}) const { - SourceVectorType source(this->source_space().mapper().size()); + VectorType source(this->source_space().mapper().size()); this->apply_inverse(range, source, opts, param); return source; } - virtual SourceVectorType - apply_inverse(const RangeVectorType& range, const std::string& type, const XT::Common::Parameter& param = {}) const + virtual VectorType + apply_inverse(const VectorType& range, const std::string& type, const XT::Common::Parameter& param = {}) const { - SourceVectorType source(this->source_space().mapper().size()); + VectorType source(this->source_space().mapper().size()); this->apply_inverse(range, source, type, param); return source; } - virtual SourceVectorType apply_inverse(const RangeVectorType& range, const XT::Common::Parameter& param = {}) const + virtual VectorType apply_inverse(const VectorType& range, const XT::Common::Parameter& param = {}) const { - SourceVectorType source(this->source_space().mapper().size()); + VectorType source(this->source_space().mapper().size()); this->apply_inverse(range, source, param); return source; } @@ -278,17 +314,55 @@ invert_options(some_type).get<std::string>("type") == some_type /// \name These jacobian variants are provided for convenience. /// \{ - virtual std::shared_ptr<ThisType> - jacobian(const SourceVectorType& source, const std::string& type, const XT::Common::Parameter& param = {}) const + virtual void jacobian(const VectorType& source, + MatrixOperatorType& jacobian_op, + const std::string& type, + const XT::Common::Parameter& param = {}) const { - return this->jacobian(source, this->jacobian_options(type), param); + return this->jacobian(source, jacobian_op, this->jacobian_options(type), param); } - virtual std::shared_ptr<ThisType> jacobian(const SourceVectorType& source, - const XT::Common::Parameter& param = {}) const - { - return this->jacobian(source, this->jacobian_options().at(0), param); - } + virtual void + jacobian(const VectorType& source, MatrixOperatorType& jacobian_op, const XT::Common::Parameter& param = {}) const + { + return this->jacobian(source, jacobian_op, this->jacobian_options().at(0), param); + } + + virtual MatrixOperatorType jacobian(const VectorType& source, + const XT::Common::Configuration& opts, + const XT::Common::Parameter& param = {}) const + { + MatrixOperatorType jacobian_op(this->source_space().grid_view(), + this->source_space(), + this->range_space(), + make_element_and_intersection_sparsity_pattern( + this->range_space(), this->source_space(), this->source_space().grid_view())); + this->jacobian(source, jacobian_op, opts, param); + return jacobian_op; + } // ... jacobian(...) + + virtual MatrixOperatorType + jacobian(const VectorType& source, const std::string& type, const XT::Common::Parameter& param = {}) const + { + MatrixOperatorType jacobian_op(this->source_space().grid_view(), + this->source_space(), + this->range_space(), + make_element_and_intersection_sparsity_pattern( + this->range_space(), this->source_space(), this->source_space().grid_view())); + this->jacobian(source, jacobian_op, type, param); + return jacobian_op; + } // ... jacobian(...) + + virtual MatrixOperatorType jacobian(const VectorType& source, const XT::Common::Parameter& param = {}) const + { + MatrixOperatorType jacobian_op(this->source_space().grid_view(), + this->source_space(), + this->range_space(), + make_element_and_intersection_sparsity_pattern( + this->range_space(), this->source_space(), this->source_space().grid_view())); + this->jacobian(source, jacobian_op, param); + return jacobian_op; + } // ... jacobian(...) /// \} /// \name These induced_norm variants are provided for convenience. @@ -298,12 +372,11 @@ invert_options(some_type).get<std::string>("type") == some_type typename = /* Only enable this method, if */ typename std::enable_if</* param is the same as XT::Common::Parameter */ ( std::is_same<ParameterType_, XT::Common::Parameter>::value) - && /* and the vector spaces defined by SourceSpaceType/SourceVectorType and */ - /* RangeSpaceType/RangeVectorType coincide. */ ( - std::is_same<SV, RV>::value&& std::is_same<SGV, RGV>::value && (s_r == r_r) - && (s_rC == r_rC) - && std::is_same<SF, RF>::value)>::type> - FieldType induced_norm(const RangeVectorType& range, const ParameterType_& param = {}) const + && /* and the vector spaces defined by SourceSpaceType/VectorType and */ + /* RangeSpaceType/VectorType coincide. */ ( + std::is_same<V, V>::value&& std::is_same<SGV, RGV>::value && (s_r == r_r) + && (s_rC == r_rC))>::type> + FieldType induced_norm(const VectorType& range, const ParameterType_& param = {}) const { using std::sqrt; return sqrt(this->apply2(range, range, param)); @@ -416,9 +489,41 @@ invert_options(some_type).get<std::string>("type") == some_type return SourceFunctionType(this->source_space(), this->apply_inverse(range.dofs().vector(), param)); } - virtual std::shared_ptr<ThisType> jacobian(const SourceFunctionType& source, - const XT::Common::Configuration& opts, - const XT::Common::Parameter& param = {}) const + virtual void jacobian(const SourceFunctionType& source, + MatrixOperatorType& jacobian_op, + const XT::Common::Configuration& opts, + const XT::Common::Parameter& param = {}) const + { + DUNE_THROW_IF(!this->source_space().contains(source), + Exceptions::operator_error, + "this->source_space() = " << this->source_space() << "\n source.space() = " << source.space()); + return this->jacobian(source.dofs().vector(), jacobian_op, opts, param); + } + + virtual void jacobian(const SourceFunctionType& source, + MatrixOperatorType& jacobian_op, + const std::string& type, + const XT::Common::Parameter& param = {}) const + { + DUNE_THROW_IF(!this->source_space().contains(source), + Exceptions::operator_error, + "this->source_space() = " << this->source_space() << "\n source.space() = " << source.space()); + return this->jacobian(source.dofs().vector(), jacobian_op, type, param); + } + + virtual void jacobian(const SourceFunctionType& source, + MatrixOperatorType& jacobian_op, + const XT::Common::Parameter& param = {}) const + { + DUNE_THROW_IF(!this->source_space().contains(source), + Exceptions::operator_error, + "this->source_space() = " << this->source_space() << "\n source.space() = " << source.space()); + return this->jacobian(source.dofs().vector(), jacobian_op, param); + } + + virtual MatrixOperatorType jacobian(const SourceFunctionType& source, + const XT::Common::Configuration& opts, + const XT::Common::Parameter& param = {}) const { DUNE_THROW_IF(!this->source_space().contains(source), Exceptions::operator_error, @@ -426,7 +531,7 @@ invert_options(some_type).get<std::string>("type") == some_type return this->jacobian(source.dofs().vector(), opts, param); } - virtual std::shared_ptr<ThisType> + virtual MatrixOperatorType jacobian(const SourceFunctionType& source, const std::string& type, const XT::Common::Parameter& param = {}) const { DUNE_THROW_IF(!this->source_space().contains(source), @@ -435,8 +540,7 @@ invert_options(some_type).get<std::string>("type") == some_type return this->jacobian(source.dofs().vector(), type, param); } - virtual std::shared_ptr<ThisType> jacobian(const SourceFunctionType& source, - const XT::Common::Parameter& param = {}) const + virtual MatrixOperatorType jacobian(const SourceFunctionType& source, const XT::Common::Parameter& param = {}) const { DUNE_THROW_IF(!this->source_space().contains(source), Exceptions::operator_error, @@ -448,11 +552,10 @@ invert_options(some_type).get<std::string>("type") == some_type typename = /* Only enable this method, if */ typename std::enable_if</* param is the same as XT::Common::Parameter */ ( std::is_same<ParameterType_, XT::Common::Parameter>::value) - && /* and the vector spaces defined by SourceSpaceType/SourceVectorType and */ - /* RangeSpaceType/RangeVectorType coincide. */ ( - std::is_same<SV, RV>::value&& std::is_same<SGV, RGV>::value && (s_r == r_r) - && (s_rC == r_rC) - && std::is_same<SF, RF>::value)>::type> + && /* and the vector spaces defined by SourceSpaceType/VectorType and */ + /* RangeSpaceType/VectorType coincide. */ ( + std::is_same<V, V>::value&& std::is_same<SGV, RGV>::value && (s_r == r_r) + && (s_rC == r_rC))>::type> FieldType induced_norm(const RangeFunctionType& range, const ParameterType_& param = {}) const { DUNE_THROW_IF(!this->range_space().contains(range), @@ -468,4 +571,6 @@ invert_options(some_type).get<std::string>("type") == some_type } // namespace GDT } // namespace Dune +#include "matrix-based.hh" + #endif // DUNE_GDT_OPERATORS_INTERFACES_HH diff --git a/dune/gdt/operators/matrix-based.hh b/dune/gdt/operators/matrix-based.hh index c013f9cae36c3e9f498d562ff5b40d0568b77f44..f94fe87fcb13571c7d9103ab3c74ab75022bca6d 100644 --- a/dune/gdt/operators/matrix-based.hh +++ b/dune/gdt/operators/matrix-based.hh @@ -22,7 +22,7 @@ #include <dune/xt/grid/type_traits.hh> #include <dune/gdt/exceptions.hh> -#include <dune/gdt/local/assembler/two-form-assemblers.hh> +#include <dune/gdt/local/assembler/bilinear-form-assemblers.hh> #include <dune/gdt/local/bilinear-forms/interfaces.hh> #include <dune/gdt/operators/interfaces.hh> #include <dune/gdt/tools/sparsity-pattern.hh> @@ -38,40 +38,22 @@ namespace GDT { * \note See OperatorInterface for a description of the template arguments. * * \sa OperatorInterface - * \sa MatrixBasedOperator + * \sa MatrixOperator */ -template <class Matrix, - class SGV, - size_t s_r = 1, - size_t s_rC = 1, - class SF = double, - class F = double, - size_t r_r = s_r, - size_t r_rC = s_rC, - class RF = double, - class RGV = SGV, - class SourceVector = typename XT::LA::Container<typename Matrix::ScalarType, Matrix::vector_type>::VectorType, - class RangeVector = SourceVector> -class ConstMatrixBasedOperator - : public OperatorInterface<SourceVector, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, RangeVector> +template <class M, class SGV, size_t s_r = 1, size_t s_rC = 1, size_t r_r = s_r, size_t r_rC = s_rC, class RGV = SGV> +class ConstMatrixOperator : public OperatorInterface<M, SGV, s_r, s_rC, r_r, r_rC, RGV> { - static_assert(XT::LA::is_matrix<Matrix>::value, ""); - - using ThisType = - ConstMatrixBasedOperator<Matrix, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, SourceVector, RangeVector>; - using BaseType = OperatorInterface<SourceVector, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, RangeVector>; + using ThisType = ConstMatrixOperator<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; + using BaseType = OperatorInterface<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; public: - using MatrixType = Matrix; - using DofFieldType = typename MatrixType::ScalarType; - + using typename BaseType::MatrixType; + using typename BaseType::VectorType; using typename BaseType::SourceSpaceType; using typename BaseType::RangeSpaceType; + using typename BaseType::MatrixOperatorType; - using typename BaseType::SourceVectorType; - using typename BaseType::RangeVectorType; - - ConstMatrixBasedOperator(const SourceSpaceType& source_spc, const RangeSpaceType& range_spc, const MatrixType& mat) + ConstMatrixOperator(const SourceSpaceType& source_spc, const RangeSpaceType& range_spc, const MatrixType& mat) : source_space_(source_spc) , range_space_(range_spc) , matrix_(mat) @@ -85,9 +67,9 @@ public: XT::Common::Exceptions::shapes_do_not_match, "matrix_.cols() = " << matrix_.cols() << "\n source_space_.mapper().size() = " << source_space_.mapper().size()); - } // ConstMatrixBasedOperator(...) + } // ConstMatrixOperator(...) - ConstMatrixBasedOperator(ThisType&& source) + ConstMatrixOperator(ThisType&& source) : source_space_(source.source_space_) , range_space_(source.range_space_) , matrix_(source.matrix_) @@ -117,9 +99,7 @@ public: using BaseType::apply; - void apply(const SourceVectorType& source, - RangeVectorType& range, - const XT::Common::Parameter& /*param*/ = {}) const override + void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& /*param*/ = {}) const override { try { matrix_.mv(source, range); @@ -146,8 +126,8 @@ public: using BaseType::apply_inverse; - void apply_inverse(const RangeVectorType& range, - SourceVectorType& source, + void apply_inverse(const VectorType& range, + VectorType& source, const XT::Common::Configuration& opts, const XT::Common::Parameter& /*param*/ = {}) const override { @@ -159,40 +139,44 @@ public: } } // ... apply_inverse(...) - std::vector<std::string> jacobian_options() const override - { - return {"direct"}; - } - - XT::Common::Configuration jacobian_options(const std::string& type) const override - { - DUNE_THROW_IF(type != jacobian_options().at(0), - Exceptions::operator_error, - "requested jacobian type is not one of the available ones!\n\n" - << "type = " - << type - << "\njacobian_options() = " - << jacobian_options()); - return {{"type", jacobian_options().at(0)}}; - } // ... jacobian_options(...) + // std::vector<std::string> jacobian_options() const override + // { + // return {"direct"}; + // } + + // XT::Common::Configuration jacobian_options(const std::string& type) const override + // { + // DUNE_THROW_IF(type != jacobian_options().at(0), + // Exceptions::operator_error, + // "requested jacobian type is not one of the available ones!\n\n" + // << "type = " + // << type + // << "\njacobian_options() = " + // << jacobian_options()); + // return {{"type", jacobian_options().at(0)}}; + // } // ... jacobian_options(...) using BaseType::jacobian; - std::shared_ptr<BaseType> jacobian(const SourceVectorType& /*source*/, - const XT::Common::Configuration& opts, - const XT::Common::Parameter& /*param*/ = {}) const override + void jacobian(const VectorType& /*source*/, + MatrixOperatorType& /*jacobian_op*/, + const XT::Common::Configuration& /*opts*/, + const XT::Common::Parameter& /*param*/ = {}) const override { - DUNE_THROW_IF( - !opts.has_key("type"), Exceptions::operator_error, "missing key 'type' in given opts!\n\nopts = " << opts); - const auto type = opts.get<std::string>("type"); - DUNE_THROW_IF(type != jacobian_options().at(0), - Exceptions::operator_error, - "requested jacobian type is not one of the available ones!\n\n" - << "type = " - << type - << "\njacobian_options() = " - << jacobian_options()); - return std::make_shared<ThisType>(source_space_, range_space_, matrix_); + // DUNE_THROW_IF( + // !opts.has_key("type"), Exceptions::operator_error, "missing key 'type' in given opts!\n\nopts = " << + // opts); + // const auto type = opts.get<std::string>("type"); + // DUNE_THROW_IF(type != jacobian_options().at(0), + // Exceptions::operator_error, + // "requested jacobian type is not one of the available ones!\n\n" + // << "type = " + // << type + // << "\njacobian_options() = " + // << jacobian_options()); + DUNE_THROW(Exceptions::operator_error, "This operator does not provide a jacobian yet!"); + // I am not sure yet how to implement this: + // if assembled, do jacobian_op.matrix() += matrix_? } // ... jacobian(...) private: @@ -200,42 +184,25 @@ private: const RangeSpaceType& range_space_; const MatrixType& matrix_; const XT::LA::Solver<MatrixType, typename SourceSpaceType::DofCommunicatorType> linear_solver_; -}; // class ConstMatrixBasedOperator - - -template <class SGV, size_t s_r, size_t s_rC, class SF, class RGV, size_t r_r, size_t r_rC, class RF, class M> -ConstMatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, - SGV, - s_r, - s_rC, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV> -make_matrix_operator(const SpaceInterface<SGV, s_r, s_rC, SF>& source_space, - const SpaceInterface<RGV, r_r, r_rC, RF>& range_space, +}; // class ConstMatrixOperator + + +template <class SGV, size_t s_r, size_t s_rC, class F, class RGV, size_t r_r, size_t r_rC, class M> +ConstMatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, SGV, s_r, s_rC, r_r, r_rC, RGV> +make_matrix_operator(const SpaceInterface<SGV, s_r, s_rC, F>& source_space, + const SpaceInterface<RGV, r_r, r_rC, F>& range_space, const XT::LA::MatrixInterface<M>& matrix) { - return ConstMatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, - SGV, - s_r, - s_rC, - SF, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>(source_space, range_space, matrix.as_imp()); -} // ... make_matrix_operator(...) + return ConstMatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, SGV, s_r, s_rC, r_r, r_rC, RGV>( + source_space, range_space, matrix.as_imp()); +} template <class GV, size_t r, size_t rC, class F, class M> -ConstMatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC, F> +ConstMatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC> make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, const XT::LA::MatrixInterface<M>& matrix) { - return ConstMatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC, F>( + return ConstMatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC>( space, space, matrix.as_imp()); } @@ -243,62 +210,44 @@ make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, const XT::LA::Ma /** * \brief Base class for linear operators which are assembled into a matrix. * - * Similar to the GlobalAssembler, we derive from the XT::Grid::Walker and povide custom append() methods to allow to - * add local element and intersection operators. In contrast to the GlobalAssembler we already hold the target matrix - * (or create one of appropriate size and given sparsity pattern), into which we want to assemble. The operator is - * assembled by walking over the given assembly_gid_view (which defaults to the one fom the given space). This allows to - * assemble an operator only on a smaller grid view than the one given from the space (similar functionality could be - * achieved by appending this operator to another walker and by providing an appropriate filter). + * We derive from the XT::Grid::Walker and povide custom append() methods to allow to add local element and intersection + * operators. In contrast to the GlobalAssembler we already hold the target matrix (or create one of appropriate size + * and given sparsity pattern), into which we want to assemble. The operator is assembled by walking over the given + * assembly_grid_view, if you want to assemble an operator only on a smaller grid view, consider to append this operator + * to another walker or provide appropriate filters. * - * \note One could achieve similar functionality by deriving from GlobalAssembler directly, which would slightly - * simplify the implementation of the append methods. However, we do not want to expose the other append methods - * of GlobalAssembler here (it should not be possible to append a local functional to an operator), but want to - * expose the ones from the XT::Grid::Walker (it should be possible to append other element or intersection - * operators or two-forms). - * - * \note See ConstMatrixBasedOperator and OperatorInterface for a description of the template arguments. + * \note See ConstMatrixOperator and OperatorInterface for a description of the template arguments. * * \sa OperatorInterface - * \sa ConstMatrixBasedOperator + * \sa ConstMatrixOperator * \sa XT::Grid::Walker - * \sa GlobalAssembler */ -template <class M, - class AssemblyGridView, - size_t s_r = 1, - size_t s_rC = 1, - class SF = double, - class SGV = AssemblyGridView, - class F = double, - size_t r_r = s_r, - size_t r_rC = s_rC, - class RF = double, - class RGV = SGV, - class SV = typename XT::LA::Container<typename M::ScalarType, M::vector_type>::VectorType, - class RV = SV> -class MatrixBasedOperator : XT::Common::StorageProvider<M>, - public ConstMatrixBasedOperator<M, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, SV, RV>, - public XT::Grid::Walker<AssemblyGridView> +template <class M, class SGV, size_t s_r = 1, size_t s_rC = 1, size_t r_r = s_r, size_t r_rC = s_rC, class RGV = SGV> +class MatrixOperator : XT::Common::StorageProvider<M>, + public ConstMatrixOperator<M, SGV, s_r, s_rC, r_r, r_rC, RGV>, + public XT::Grid::Walker<SGV> { - static_assert(XT::Grid::is_view<AssemblyGridView>::value, ""); - - using ThisType = MatrixBasedOperator<M, AssemblyGridView, s_r, s_rC, SF, SGV, F, r_r, r_rC, RF, RGV, SV, RV>; + using ThisType = MatrixOperator<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; using MatrixStorage = XT::Common::StorageProvider<M>; - using OperatorBaseType = ConstMatrixBasedOperator<M, SGV, s_r, s_rC, SF, F, r_r, r_rC, RF, RGV, SV, RV>; - using WalkerBaseType = XT::Grid::Walker<AssemblyGridView>; + using OperatorBaseType = ConstMatrixOperator<M, SGV, s_r, s_rC, r_r, r_rC, RGV>; + using WalkerBaseType = XT::Grid::Walker<SGV>; public: - using AssemblyGridViewType = AssemblyGridView; + using AssemblyGridViewType = SGV; using typename OperatorBaseType::MatrixType; + using typename OperatorBaseType::VectorType; using typename OperatorBaseType::SourceSpaceType; using typename OperatorBaseType::RangeSpaceType; + using typename OperatorBaseType::MatrixOperatorType; + using typename OperatorBaseType::V; + using typename OperatorBaseType::F; using typename WalkerBaseType::ElementType; - using ElementFilterType = XT::Grid::ElementFilter<AssemblyGridViewType>; - using IntersectionFilterType = XT::Grid::IntersectionFilter<AssemblyGridViewType>; - using ApplyOnAllElements = XT::Grid::ApplyOn::AllElements<AssemblyGridViewType>; - using ApplyOnAllIntersections = XT::Grid::ApplyOn::AllIntersections<AssemblyGridViewType>; + using ElementFilterType = XT::Grid::ElementFilter<SGV>; + using IntersectionFilterType = XT::Grid::IntersectionFilter<SGV>; + using ApplyOnAllElements = XT::Grid::ApplyOn::AllElements<SGV>; + using ApplyOnAllIntersections = XT::Grid::ApplyOn::AllIntersections<SGV>; using typename WalkerBaseType::E; using typename WalkerBaseType::I; @@ -306,10 +255,10 @@ public: /** * Ctor which accept an existing matrix into which to assemble. */ - MatrixBasedOperator(AssemblyGridViewType assembly_grid_view, - const SourceSpaceType& source_spc, - const RangeSpaceType& range_spc, - MatrixType& mat) + MatrixOperator(AssemblyGridViewType assembly_grid_view, + const SourceSpaceType& source_spc, + const RangeSpaceType& range_spc, + MatrixType& mat) : MatrixStorage(mat) , OperatorBaseType(source_spc, range_spc, MatrixStorage::access()) , WalkerBaseType(assembly_grid_view) @@ -324,10 +273,10 @@ public: * Ctor which creates an appropriate matrix into which to assemble from a given sparsity pattern. */ - MatrixBasedOperator(AssemblyGridViewType assembly_grid_view, - const SourceSpaceType& source_spc, - const RangeSpaceType& range_spc, - const XT::LA::SparsityPatternDefault& pattern) + MatrixOperator(AssemblyGridViewType assembly_grid_view, + const SourceSpaceType& source_spc, + const RangeSpaceType& range_spc, + const XT::LA::SparsityPatternDefault& pattern) : MatrixStorage(new MatrixType(range_spc.mapper().size(), source_spc.mapper().size(), pattern)) , OperatorBaseType(source_spc, range_spc, MatrixStorage::access()) , WalkerBaseType(assembly_grid_view) @@ -347,40 +296,29 @@ public: using WalkerBaseType::append; - ThisType& append(const LocalElementBilinearFormInterface<E, r_r, r_rC, RF, F, s_r, s_rC, SF>& local_bilinear_form, + ThisType& append(const LocalElementBilinearFormInterface<E, r_r, r_rC, F, F, s_r, s_rC, F>& local_bilinear_form, const XT::Common::Parameter& param = {}, const ElementFilterType& filter = ApplyOnAllElements()) { - using LocalAssemblerType = - LocalElementBilinearFormAssembler<MatrixType, AssemblyGridViewType, r_r, r_rC, RF, RGV, SGV, s_r, s_rC, SF>; + using LocalAssemblerType = LocalElementBilinearFormAssembler<M, SGV, r_r, r_rC, F, RGV, SGV, s_r, s_rC, F>; this->append(new LocalAssemblerType( this->range_space(), this->source_space(), local_bilinear_form, MatrixStorage::access(), param), filter); return *this; } - ThisType& - append(const LocalIntersectionBilinearFormInterface<I, r_r, r_rC, RF, F, s_r, s_rC, SF>& local_bilinear_form, - const XT::Common::Parameter& param = {}, - const IntersectionFilterType& filter = ApplyOnAllIntersections()) + ThisType& append(const LocalIntersectionBilinearFormInterface<I, r_r, r_rC, F, F, s_r, s_rC, F>& local_bilinear_form, + const XT::Common::Parameter& param = {}, + const IntersectionFilterType& filter = ApplyOnAllIntersections()) { - using LocalAssemblerType = LocalIntersectionBilinearFormAssembler<MatrixType, - AssemblyGridViewType, - r_r, - r_rC, - RF, - RGV, - SGV, - s_r, - s_rC, - SF>; + using LocalAssemblerType = + LocalIntersectionBilinearFormAssembler<MatrixType, AssemblyGridViewType, r_r, r_rC, F, RGV, SGV, s_r, s_rC, F>; this->append(new LocalAssemblerType( this->range_space(), this->source_space(), local_bilinear_form, MatrixStorage::access(), param), filter); return *this; } // ... append(...) - // similar append for LocalIntersectionFunctionalInterface ... void assemble(const bool use_tbb = false) override final { @@ -391,68 +329,52 @@ public: assembled_ = true; } + using OperatorBaseType::jacobian; + + void jacobian(const VectorType& /*source*/, + MatrixOperatorType& /*jacobian_op*/, + const XT::Common::Configuration& /*opts*/, + const XT::Common::Parameter& /*param*/ = {}) const override + { + DUNE_THROW(Exceptions::operator_error, "This operator does not provide a jacobian yet!"); + // I am not sure yet how to implement this: + // if assembled, do jacobian_op.matrix() += matrix_? + // If not, append this to jacobian_op? + // Or keep a list of all appended local ops and append them to jacobian_op? <- This is probably best. + } // ... jacobian(...) + private: bool assembled_; -}; // class MatrixBasedOperator +}; // class MatrixOperator /// \name Variants of make_matrix_operator for a given matrix. /// \{ -template <class AGV, - class SGV, - size_t s_r, - size_t s_rC, - class SF, - class RGV, - size_t r_r, - size_t r_rC, - class RF, - class M> -MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV> -make_matrix_operator(GridView<AGV> assembly_grid_view, - const SpaceInterface<SGV, s_r, s_rC, SF>& source_space, - const SpaceInterface<RGV, r_r, r_rC, RF>& range_space, +template <class SGV, size_t s_r, size_t s_rC, class F, class RGV, size_t r_r, size_t r_rC, class M> +MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, SGV, s_r, s_rC, r_r, r_rC, RGV> +make_matrix_operator(SGV assembly_grid_view, + const SpaceInterface<SGV, s_r, s_rC, F>& source_space, + const SpaceInterface<RGV, r_r, r_rC, F>& range_space, XT::LA::MatrixInterface<M>& matrix) { - return MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>(assembly_grid_view, source_space, range_space, matrix.as_imp()); -} // ... make_matrix_operator(...) - -template <class AGV, class GV, size_t r, size_t rC, class F, class M> -MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GridView<AGV>, r, rC, F, GV> -make_matrix_operator(GridView<AGV> assembly_grid_view, - const SpaceInterface<GV, r, rC, F>& space, - XT::LA::MatrixInterface<M>& matrix) + return MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, SGV, s_r, s_rC, r_r, r_rC, RGV>( + assembly_grid_view, source_space, range_space, matrix.as_imp()); +} + +template <class GV, size_t r, size_t rC, class F, class M> +MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC> make_matrix_operator( + GV assembly_grid_view, const SpaceInterface<GV, r, rC, F>& space, XT::LA::MatrixInterface<M>& matrix) { - return MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GridView<AGV>, r, rC, F, GV>( + return MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC>( assembly_grid_view, space, space, matrix.as_imp()); } template <class GV, size_t r, size_t rC, class F, class M> -MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC, F> +MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC> make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, XT::LA::MatrixInterface<M>& matrix) { - return MatrixBasedOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC, F>( + return MatrixOperator<typename XT::LA::MatrixInterface<M>::derived_type, GV, r, rC>( space.grid_view(), space, space, matrix.as_imp()); } @@ -466,45 +388,17 @@ make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, XT::LA::MatrixIn auto op = make_matrix_operator<MatrixType>(assembly_grid_view, source_space, range_space, pattern); \endcode */ -template <class MatrixType, - class AGV, - class SGV, - size_t s_r, - size_t s_rC, - class SF, - class RGV, - size_t r_r, - size_t r_rC, - class RF> +template <class MatrixType, class SGV, size_t s_r, size_t s_rC, class F, class RGV, size_t r_r, size_t r_rC> typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, - MatrixBasedOperator<MatrixType, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>>::type -make_matrix_operator(GridView<AGV> assembly_grid_view, - const SpaceInterface<SGV, s_r, s_rC, SF>& source_space, - const SpaceInterface<RGV, r_r, r_rC, RF>& range_space, + MatrixOperator<MatrixType, SGV, s_r, s_rC, r_r, r_rC, RGV>>::type +make_matrix_operator(SGV assembly_grid_view, + const SpaceInterface<SGV, s_r, s_rC, F>& source_space, + const SpaceInterface<RGV, r_r, r_rC, F>& range_space, const XT::LA::SparsityPatternDefault& pattern) { - return MatrixBasedOperator<MatrixType, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>(assembly_grid_view, source_space, range_space, pattern); -} // ... make_matrix_operator(...) + return MatrixOperator<MatrixType, SGV, s_r, s_rC, r_r, r_rC, RGV>( + assembly_grid_view, source_space, range_space, pattern); +} /** * \note Use as in @@ -512,14 +406,13 @@ make_matrix_operator(GridView<AGV> assembly_grid_view, auto op = make_matrix_operator<MatrixType>(assembly_grid_view, space, pattern); \endcode */ -template <class MatrixType, class AGV, class GV, size_t r, size_t rC, class F> -typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, - MatrixBasedOperator<MatrixType, GridView<AGV>, r, rC, F, GV>>::type -make_matrix_operator(GridView<AGV> assembly_grid_view, +template <class MatrixType, class GV, size_t r, size_t rC, class F> +typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixOperator<MatrixType, GV, r, rC>>::type +make_matrix_operator(GV assembly_grid_view, const SpaceInterface<GV, r, rC, F>& space, const XT::LA::SparsityPatternDefault& pattern) { - return MatrixBasedOperator<MatrixType, GridView<AGV>, r, rC, F, GV>(assembly_grid_view, space, space, pattern); + return MatrixOperator<MatrixType, GV, r, rC>(assembly_grid_view, space, space, pattern); } /** @@ -529,10 +422,10 @@ auto op = make_matrix_operator<MatrixType>(space, pattern); \endcode */ template <class MatrixType, class GV, size_t r, size_t rC, class F> -typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixBasedOperator<MatrixType, GV, r, rC, F>>::type +typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixOperator<MatrixType, GV, r, rC>>::type make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, const XT::LA::SparsityPatternDefault& pattern) { - return MatrixBasedOperator<MatrixType, GV, r, rC, F>(space.grid_view(), space, space, pattern); + return MatrixOperator<MatrixType, GV, r, rC>(space.grid_view(), space, space, pattern); } /// \} @@ -545,48 +438,20 @@ make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, const XT::LA::Sp auto op = make_matrix_operator<MatrixType>(source_space, range_space, stencil); \endcode */ -template <class MatrixType, - class AGV, - class SGV, - size_t s_r, - size_t s_rC, - class SF, - class RGV, - size_t r_r, - size_t r_rC, - class RF> +template <class MatrixType, class SGV, size_t s_r, size_t s_rC, class F, class RGV, size_t r_r, size_t r_rC> typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, - MatrixBasedOperator<MatrixType, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>>::type -make_matrix_operator(GridView<AGV> assembly_grid_view, - const SpaceInterface<SGV, s_r, s_rC, SF>& source_space, - const SpaceInterface<RGV, r_r, r_rC, RF>& range_space, + MatrixOperator<MatrixType, SGV, s_r, s_rC, r_r, r_rC, RGV>>::type +make_matrix_operator(SGV assembly_grid_view, + const SpaceInterface<SGV, s_r, s_rC, F>& source_space, + const SpaceInterface<RGV, r_r, r_rC, F>& range_space, const Stencil stencil = Stencil::element_and_intersection) { - return MatrixBasedOperator<MatrixType, - GridView<AGV>, - s_r, - s_rC, - SF, - SGV, - typename XT::Common::multiplication_promotion<SF, RF>::type, - r_r, - r_rC, - RF, - RGV>(assembly_grid_view, - source_space, - range_space, - make_sparsity_pattern(range_space, source_space, assembly_grid_view, stencil)); -} // ... make_matrix_operator(...) + return MatrixOperator<MatrixType, SGV, s_r, s_rC, r_r, r_rC, RGV>( + assembly_grid_view, + source_space, + range_space, + make_sparsity_pattern(range_space, source_space, assembly_grid_view, stencil)); +} /** * \note Use as in @@ -594,14 +459,13 @@ make_matrix_operator(GridView<AGV> assembly_grid_view, auto op = make_matrix_operator<MatrixType>(assembly_grid_view, space, stencil); \endcode */ -template <class MatrixType, class AGV, class GV, size_t r, size_t rC, class F> -typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, - MatrixBasedOperator<MatrixType, GridView<AGV>, r, rC, F, GV>>::type -make_matrix_operator(GridView<AGV> assembly_grid_view, +template <class MatrixType, class GV, size_t r, size_t rC, class F> +typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixOperator<MatrixType, GV, r, rC>>::type +make_matrix_operator(GV assembly_grid_view, const SpaceInterface<GV, r, rC, F>& space, const Stencil stencil = Stencil::element_and_intersection) { - return MatrixBasedOperator<MatrixType, GridView<AGV>, r, rC, F, GV>( + return MatrixOperator<MatrixType, GV, r, rC>( assembly_grid_view, space, space, make_sparsity_pattern(space, assembly_grid_view, stencil)); } @@ -612,12 +476,11 @@ auto op = make_matrix_operator<MatrixType>(space, stencil); \endcode */ template <class MatrixType, class GV, size_t r, size_t rC, class F> -typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixBasedOperator<MatrixType, GV, r, rC, F>>::type +typename std::enable_if<XT::LA::is_matrix<MatrixType>::value, MatrixOperator<MatrixType, GV, r, rC>>::type make_matrix_operator(const SpaceInterface<GV, r, rC, F>& space, const Stencil stencil = Stencil::element_and_intersection) { - return MatrixBasedOperator<MatrixType, GV, r, rC, F>( - space.grid_view(), space, space, make_sparsity_pattern(space, stencil)); + return MatrixOperator<MatrixType, GV, r, rC>(space.grid_view(), space, space, make_sparsity_pattern(space, stencil)); } /// \}