From 86de7d1603e3e46a2a7bd581a7d9ee646872bea1 Mon Sep 17 00:00:00 2001
From: Felix Schindler <felix.schindler@wwu.de>
Date: Tue, 7 Aug 2018 15:19:07 +0200
Subject: [PATCH] [operators.advection-fv] add boundary treatment by custom
 numerical flux

---
 dune/gdt/local/operators/advection-fv.hh | 24 +++++++------
 dune/gdt/operators/advection-fv.hh       | 44 ++++++++++++++++++++----
 2 files changed, 50 insertions(+), 18 deletions(-)

diff --git a/dune/gdt/local/operators/advection-fv.hh b/dune/gdt/local/operators/advection-fv.hh
index 77d2fd15d..d1531c64a 100644
--- a/dune/gdt/local/operators/advection-fv.hh
+++ b/dune/gdt/local/operators/advection-fv.hh
@@ -532,9 +532,11 @@ public:
     const auto normal = intersection.centerUnitOuterNormal();
     const auto g = numerical_flux_->apply(u->dofs(), v->dofs(), normal, param);
     const auto h_intersection = intersection.geometry().volume();
+    const auto h_inside_element = inside_element.geometry().volume();
+    const auto h_outside_element = outside_element.geometry().volume();
     for (size_t ii = 0; ii < m; ++ii) {
-      local_range_inside.dofs()[ii] += (g[ii] * h_intersection) / inside_element.geometry().volume();
-      local_range_outside.dofs()[ii] -= (g[ii] * h_intersection) / outside_element.geometry().volume();
+      local_range_inside.dofs()[ii] += (g[ii] * h_intersection) / h_inside_element;
+      local_range_outside.dofs()[ii] -= (g[ii] * h_intersection) / h_outside_element;
     }
   } // ... apply(...)
 
@@ -544,10 +546,10 @@ private:
 
 
 template <class I, class SV, class SGV, size_t m = 1, class SF = double, class RF = SF, class RGV = SGV, class RV = SV>
-class LocalAdvectionFvBoundaryOperatorByCustomNumericalFlux
+class LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator
     : public LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SF, m, 1, RF, RGV, RV>
 {
-  using ThisType = LocalAdvectionFvBoundaryOperatorByCustomNumericalFlux<I, SV, SGV, m, SF, RF, RGV, RV>;
+  using ThisType = LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator<I, SV, SGV, m, SF, RF, RGV, RV>;
   using BaseType = LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SF, m, 1, RF, RGV, RV>;
 
 public:
@@ -563,16 +565,16 @@ public:
   using LambdaType = std::function<StateRangeType(
       const StateDofsType& /*u*/, const StateDomainType& /*n*/, const XT::Common::Parameter& /*param*/)>;
 
-  LocalAdvectionFvBoundaryOperatorByCustomNumericalFlux(
+  LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
       LambdaType numerical_boundary_flux_lambda, const XT::Common::ParameterType& boundary_treatment_param_type = {})
     : BaseType(boundary_treatment_param_type)
-    , numerical_boundary_flux_lambda_(numerical_boundary_flux_lambda)
+    , numerical_boundary_flux_(numerical_boundary_flux_lambda)
   {
   }
 
-  LocalAdvectionFvBoundaryOperatorByCustomNumericalFlux(const ThisType& other)
+  LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(const ThisType& other)
     : BaseType(other.parameter_type())
-    , numerical_boundary_flux_lambda_(other.numerical_boundary_flux_lambda_)
+    , numerical_boundary_flux_(other.numerical_boundary_flux_)
   {
   }
 
@@ -594,7 +596,7 @@ public:
     const auto& element = local_range_inside.element();
     const auto u = source.local_discrete_function(element);
     const auto normal = intersection.centerUnitOuterNormal();
-    const auto g = numerical_boundary_flux_lambda_(u->dofs(), normal, param);
+    const auto g = numerical_boundary_flux_(u->dofs(), normal, param);
     const auto h_intersection = intersection.geometry().volume();
     const auto h_element = element.geometry().volume();
     for (size_t ii = 0; ii < m; ++ii)
@@ -602,8 +604,8 @@ public:
   } // ... apply(...)
 
 private:
-  const LambdaType numerical_boundary_flux_lambda_;
-}; // class LocalAdvectionFvBoundaryOperatorByCustomNumericalFlux
+  const LambdaType numerical_boundary_flux_;
+}; // class LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator
 
 
 } // namespace GDT
diff --git a/dune/gdt/operators/advection-fv.hh b/dune/gdt/operators/advection-fv.hh
index b4d0f077a..2859eba46 100644
--- a/dune/gdt/operators/advection-fv.hh
+++ b/dune/gdt/operators/advection-fv.hh
@@ -12,6 +12,7 @@
 
 #include <dune/xt/common/type_traits.hh>
 #include <dune/xt/grid/type_traits.hh>
+#include <dune/xt/grid/walker/filters.hh>
 
 #include <dune/gdt/local/operators/advection-fv.hh>
 
@@ -69,18 +70,22 @@ public:
   using AssemblyGridViewType = AssemblyGridView;
   using NumericalFluxType = NumericalFluxInterface<AssemblyGridView::dimension, m, RF>;
 
+  using I = XT::Grid::extract_intersection_t<AGV>;
+  using BoundaryTreatmentByCustomNumericalFluxOperatorType =
+      LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator<I, SV, SGV, m, SF, RF, RGV, RV>;
+
   using typename BaseType::SourceSpaceType;
   using typename BaseType::RangeSpaceType;
 
   using typename BaseType::SourceVectorType;
   using typename BaseType::RangeVectorType;
 
-  AdvectionFvOperator(const AssemblyGridViewType& assembly_grid_view,
-                      const NumericalFluxType& numerical_flux,
-                      const SourceSpaceType& source_space,
-                      const RangeSpaceType& range_space,
-                      const XT::Grid::IntersectionFilter<AssemblyGridViewType>& periodicity_exception =
-                          XT::Grid::ApplyOn::NoIntersections<AssemblyGridViewType>())
+  AdvectionFvOperator(
+      const AssemblyGridViewType& 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>())
     : BaseType(numerical_flux.parameter_type())
     , assembly_grid_view_(assembly_grid_view)
     , numerical_flux_(numerical_flux)
@@ -107,6 +112,23 @@ public:
     return range_space_;
   }
 
+  /// \name Non-periodic boundary treatment
+  /// \{
+
+  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>())
+  {
+    boundary_treatments_by_custom_numerical_flux_.emplace_back(
+        new BoundaryTreatmentByCustomNumericalFluxOperatorType(numerical_boundary_treatment_flux,
+                                                               boundary_treatment_parameter_type),
+        filter.copy());
+    return *this;
+  }
+
+  /// \}
+
   using BaseType::apply;
 
   void apply(const SourceVectorType& source,
@@ -123,7 +145,6 @@ public:
     auto range_function = make_discrete_function(range_space_, range);
     // set up the actual operator
     auto localizable_op = make_localizable_operator(assembly_grid_view_, source_function, range_function);
-    using I = XT::Grid::extract_intersection_t<AGV>;
     // treat all inner intersections
     localizable_op.append(LocalAdvectionFvCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(numerical_flux_),
                           param,
@@ -132,6 +153,12 @@ public:
     localizable_op.append(LocalAdvectionFvCouplingOperator<I, SV, SGV, m, SF, RF, RGV, RV>(numerical_flux_),
                           param,
                           *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<AGV>() && !(*periodicity_exception_)));
+    // treat boundaries by custom numerical flux
+    for (const auto& boundary_treatment : boundary_treatments_by_custom_numerical_flux_) {
+      const auto& boundary_op = *boundary_treatment.first;
+      const auto& filter = *boundary_treatment.second;
+      localizable_op.append(boundary_op, param, filter);
+    }
     // do the actual work
     localizable_op.assemble(/*use_tbb=*/true);
     DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!");
@@ -143,6 +170,9 @@ private:
   const SourceSpaceType& source_space_;
   const RangeSpaceType& range_space_;
   std::unique_ptr<XT::Grid::IntersectionFilter<AssemblyGridViewType>> periodicity_exception_;
+  std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>,
+                      std::unique_ptr<XT::Grid::IntersectionFilter<AGV>>>>
+      boundary_treatments_by_custom_numerical_flux_;
 }; // class AdvectionFvOperator
 
 
-- 
GitLab