diff --git a/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh b/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
index c5929c639cfb08838616a6516f62b6eeeb35f104..64082860486db731a1e26c990f7402067280c298 100644
--- a/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
@@ -11,10 +11,11 @@
#ifndef DUNE_GDT_HYPERBOLIC_PROBLEMS_MOMENTMODELS_HATFUNCTIONS_HH
#define DUNE_GDT_HYPERBOLIC_PROBLEMS_MOMENTMODELS_HATFUNCTIONS_HH
-#include <vector>
#include <string>
+#include <vector>
#include <dune/xt/common/fmatrix.hh>
+#include <dune/xt/common/numeric.hh>
#include <dune/gdt/test/momentmodels/triangulation.hh>
@@ -24,6 +25,129 @@ namespace Dune {
namespace GDT {
+template <class DomainFieldType,
+ size_t domainDim,
+ class RangeFieldType,
+ size_t rangeDim,
+ size_t fluxDim,
+ EntropyType entropy>
+class HatFunctionMomentBasisBase
+ : public MomentBasisInterface<DomainFieldType, domainDim, RangeFieldType, rangeDim, 1, fluxDim, entropy>
+{
+ using BaseType = MomentBasisInterface<DomainFieldType, domainDim, RangeFieldType, rangeDim, 1, fluxDim, entropy>;
+
+public:
+ using BaseType::dimDomain;
+ using BaseType::dimRange;
+
+ using typename BaseType::DynamicRangeType;
+ using typename BaseType::RangeType;
+ using typename BaseType::StringifierType;
+
+ template <class... Args>
+ HatFunctionMomentBasisBase(Args&&... args)
+ : BaseType(std::forward<Args>(args)...)
+ {}
+
+ static StringifierType stringifier()
+ {
+ return [](const RangeType& val) {
+ RangeFieldType psi(0);
+ for (const auto& entry : val)
+ psi += entry;
+ return XT::Common::to_string(psi, 15);
+ };
+ } // ... stringifier()
+
+ DynamicRangeType alpha_one() const override final
+ {
+ return DynamicRangeType(dimRange, 1.);
+ }
+
+ RangeFieldType density(const DynamicRangeType& u) const override final
+ {
+ return XT::Common::reduce(&(u[0]), &(u[0]) + dimRange, RangeFieldType(0));
+ }
+
+ RangeFieldType density(const RangeType& u) const override final
+ {
+ return XT::Common::reduce(&(u[0]), &(u[0]) + dimRange, RangeFieldType(0));
+ }
+
+ template <class Vec>
+ std::enable_if_t<XT::Common::is_vector<Vec>::value, void> alpha_one(Vec& ret) const
+ {
+ for (size_t ii = 0; ii < ret.size(); ++ii)
+ XT::Common::VectorAbstraction<Vec>::set_entry(ret, ii, 1.);
+ }
+
+ template <class Vec>
+ std::enable_if_t<XT::Common::is_vector<Vec>::value && !std::is_same<Vec, DynamicRangeType>::value, RangeFieldType>
+ density(const Vec& u) const
+ {
+ RangeFieldType ret(0.);
+ for (size_t ii = 0; ii < u.size(); ++ii)
+ ret += XT::Common::VectorAbstraction<Vec>::get_entry(u, ii);
+ return ret;
+ }
+
+ using BaseType::u_iso;
+
+ template <class Vec>
+ std::enable_if_t<XT::Common::is_vector<Vec>::value, void> u_iso(Vec& ret) const
+ {
+ const auto& ret_range = u_iso();
+ using V = XT::Common::VectorAbstraction<Vec>;
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ V::set_entry(ret, ii, ret_range[ii]);
+ }
+
+ void ensure_min_density(DynamicRangeType& u, const RangeFieldType psi_min) const override final
+ {
+ const auto u_iso_min = u_iso() * psi_min;
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ if (u[ii] < u_iso_min[ii])
+ u[ii] = u_iso_min[ii];
+ }
+
+ void ensure_min_density(RangeType& u, const RangeFieldType psi_min) const override final
+ {
+ const auto u_iso_min = u_iso() * psi_min;
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ if (u[ii] < u_iso_min[ii])
+ u[ii] = u_iso_min[ii];
+ }
+
+ bool adjust_alpha_to_ensure_min_density(RangeType& alpha, const RangeFieldType psi_min) const override final
+ {
+ bool changed = false;
+ const auto alpha_min = std::log(psi_min);
+ for (size_t ii = 0; ii < dimRange; ++ii) {
+ if (alpha[ii] < alpha_min) {
+ alpha[ii] = alpha_min;
+ changed = true;
+ }
+ }
+ return changed;
+ }
+
+ std::string short_id() const override final
+ {
+ return "hf";
+ }
+
+ std::string mn_name() const override final
+ {
+ return "hfm" + XT::Common::to_string(dimRange);
+ }
+
+ std::string pn_name() const override final
+ {
+ return "hfp" + XT::Common::to_string(dimRange);
+ }
+}; // class HatFunctionMomentBasisBase
+
+
template <class DomainFieldType,
size_t dimDomain,
class RangeFieldType,
@@ -43,7 +167,7 @@ template <class DomainFieldType,
size_t fluxDim,
EntropyType entropy>
class HatFunctionMomentBasis<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, fluxDim, entropy>
- : public MomentBasisInterface<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, fluxDim, entropy>
+ : public HatFunctionMomentBasisBase<DomainFieldType, 1, RangeFieldType, rangeDim, fluxDim, entropy>
{
public:
static const size_t dimDomain = 1;
@@ -52,8 +176,7 @@ public:
static const size_t num_intervals = dimRange - 1;
private:
- using BaseType =
- MomentBasisInterface<DomainFieldType, dimDomain, RangeFieldType, dimRange, dimRangeCols, fluxDim, entropy>;
+ using BaseType = HatFunctionMomentBasisBase<DomainFieldType, 1, RangeFieldType, rangeDim, fluxDim, entropy>;
public:
using typename BaseType::DomainType;
@@ -63,14 +186,7 @@ public:
using typename BaseType::RangeType;
using typename BaseType::SphericalTriangulationType;
using typename BaseType::StringifierType;
- using typename BaseType::VisualizerType;
using PartitioningType = typename BaseType::Partitioning1dType;
- template <class DiscreteFunctionType>
-
- static std::string static_id()
- {
- return "hatfunctions";
- }
static size_t default_quad_order()
{
@@ -130,18 +246,6 @@ public:
return ret;
} // ... evaluate(...)
- DynamicRangeType integrated() const override final
- {
- DynamicRangeType ret(dimRange, 0);
- const auto& mu = partitioning_;
- ret[0] = mu[1] - mu[0];
- for (size_t ii = 1; ii < num_intervals; ++ii)
- ret[ii] = mu[ii + 1] - mu[ii - 1];
- ret[num_intervals] = mu[num_intervals] - mu[num_intervals - 1];
- ret *= 0.5;
- return ret;
- }
-
// returns matrix with entries <h_i h_j>
MatrixType mass_matrix() const override final
{
@@ -259,76 +363,58 @@ public:
return ret;
}
- static StringifierType stringifier()
- {
- return [](const RangeType& val) {
- RangeFieldType psi(0);
- for (const auto& entry : val)
- psi += entry;
- return XT::Common::to_string(psi, 15);
- };
- } // ... stringifier()
-
- const PartitioningType& partitioning() const
- {
- return partitioning_;
- }
-
- DynamicRangeType alpha_one() const override final
- {
- return DynamicRangeType(dimRange, 1.);
- }
-
- RangeFieldType density(const DynamicRangeType& u) const override final
- {
- return std::accumulate(u.begin(), u.end(), RangeFieldType(0));
- }
-
- using BaseType::u_iso;
-
- void ensure_min_density(DynamicRangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t ii = 0; ii < dimRange; ++ii)
- if (u[ii] < u_iso_min[ii])
- u[ii] = u_iso_min[ii];
- }
-
- void ensure_min_density(RangeType& u, const RangeFieldType min_density) const override final
+ // get indices of all faces that contain point v
+ std::vector<size_t> get_face_indices(const DomainType& v) const
{
- const auto u_iso_min = u_iso() * min_density;
- for (size_t ii = 0; ii < dimRange; ++ii)
- if (u[ii] < u_iso_min[ii])
- u[ii] = u_iso_min[ii];
+ std::vector<size_t> face_indices;
+ for (size_t jj = 0; jj < partitioning_.size() - 1; ++jj)
+ if (XT::Common::FloatCmp::ge(v[0], partitioning_[jj]) && XT::Common::FloatCmp::le(v[0], partitioning_[jj + 1]))
+ face_indices.push_back(jj);
+ assert(face_indices.size());
+ return face_indices;
}
- std::string short_id() const override final
+ DynamicRangeType integrated_initializer(const QuadraturesType& /*quadratures*/) const override final
{
- return "hf";
+ DynamicRangeType ret(dimRange, 0);
+ const auto& mu = partitioning_;
+ ret[0] = mu[1] - mu[0];
+ for (size_t ii = 1; ii < num_intervals; ++ii)
+ ret[ii] = mu[ii + 1] - mu[ii - 1];
+ ret[num_intervals] = mu[num_intervals] - mu[num_intervals - 1];
+ ret *= 0.5;
+ return ret;
}
- std::string mn_name() const override final
+ const PartitioningType& partitioning() const
{
- return "hfm" + XT::Common::to_string(dimRange);
+ return partitioning_;
}
- std::string pn_name() const override final
+ DynamicRangeType
+ get_u(const std::function<RangeFieldType(DomainType, std::array<int, dimDomain>)>& psi) const override final
{
- return "hfp" + XT::Common::to_string(dimRange);
+ DynamicRangeType ret(dimRange, 0.);
+ const auto merged_quads = XT::Data::merged_quadrature(quadratures_);
+ for (auto it = merged_quads.begin(); it != merged_quads.end(); ++it) {
+ const auto& quad_point = *it;
+ const size_t interval_index = it.first_index();
+ const auto& v = quad_point.position();
+ const auto val = evaluate_on_interval(v, interval_index);
+ const auto factor = psi(v, has_fixed_sign(interval_index)) * quad_point.weight();
+ for (size_t ii = 0; ii < 2; ++ii)
+ ret[interval_index + ii] += val[ii] * factor;
+ }
+ return ret;
}
- // get indices of all faces that contain point v
- std::vector<size_t> get_face_indices(const DomainType& v) const
+ std::array<int, dimDomain> has_fixed_sign(const size_t index) const override final
{
- std::vector<size_t> face_indices;
- for (size_t jj = 0; jj < partitioning_.size() - 1; ++jj)
- if (XT::Common::FloatCmp::ge(v[0], partitioning_[jj]) && XT::Common::FloatCmp::le(v[0], partitioning_[jj + 1]))
- face_indices.push_back(jj);
- assert(face_indices.size());
- return face_indices;
+ return BaseType::interval_has_fixed_sign(index, num_intervals);
}
private:
+ using BaseType::quadratures_;
const PartitioningType partitioning_;
}; // class HatFunctionMomentBasis<DomainFieldType, 1, ...>
@@ -343,13 +429,12 @@ constexpr size_t
template <class DomainFieldType, class RangeFieldType, size_t refinements, size_t fluxDim, EntropyType entropy>
class HatFunctionMomentBasis<DomainFieldType, 3, RangeFieldType, refinements, 1, fluxDim, entropy>
- : public MomentBasisInterface<DomainFieldType,
- 3,
- RangeFieldType,
- OctaederStatistics<refinements>::num_vertices(),
- 1,
- fluxDim,
- entropy>
+ : public HatFunctionMomentBasisBase<DomainFieldType,
+ 3,
+ RangeFieldType,
+ OctaederStatistics<refinements>::num_vertices(),
+ fluxDim,
+ entropy>
{
public:
static constexpr size_t dimDomain = 3;
@@ -359,8 +444,12 @@ public:
static constexpr size_t num_refinements = refinements;
private:
- using BaseType =
- MomentBasisInterface<DomainFieldType, dimDomain, RangeFieldType, dimRange, dimRangeCols, dimFlux, entropy>;
+ using BaseType = HatFunctionMomentBasisBase<DomainFieldType,
+ 3,
+ RangeFieldType,
+ OctaederStatistics<refinements>::num_vertices(),
+ fluxDim,
+ entropy>;
using ThisType = HatFunctionMomentBasis;
public:
@@ -372,11 +461,9 @@ public:
using typename BaseType::QuadraturesType;
using typename BaseType::RangeType;
using typename BaseType::StringifierType;
- using typename BaseType::VisualizerType;
using LocalMatrixType = XT::Common::FieldMatrix<RangeFieldType, 3, 3>;
using BaseType::barycentre_rule;
- using BaseType::is_negative;
static size_t default_quad_order()
{
@@ -461,16 +548,6 @@ public:
return ret;
} // ... evaluate(...)
- static StringifierType stringifier()
- {
- return [](const RangeType& val) {
- RangeFieldType psi(0);
- for (const auto& entry : val)
- psi += entry;
- return XT::Common::to_string(psi, 15);
- };
- } // ... stringifier()
-
const TriangulationType& triangulation() const
{
return triangulation_;
@@ -481,48 +558,6 @@ public:
return BaseType::unit_ball_volume_quad();
}
- DynamicRangeType alpha_one() const override final
- {
- return DynamicRangeType(dimRange, 1.);
- }
-
- template <class Vec>
- std::enable_if_t<XT::Common::is_vector<Vec>::value, void> alpha_one(Vec& ret) const
- {
- for (size_t ii = 0; ii < ret.size(); ++ii)
- XT::Common::VectorAbstraction<Vec>::set_entry(ret, ii, 1.);
- }
-
- RangeFieldType density(const DynamicRangeType& u) const override final
- {
- return std::accumulate(u.begin(), u.end(), 0.);
- }
-
- template <class Vec>
- std::enable_if_t<XT::Common::is_vector<Vec>::value && !std::is_same<Vec, DynamicRangeType>::value, RangeFieldType>
- density(const Vec& u) const
- {
- RangeFieldType ret(0.);
- for (size_t ii = 0; ii < u.size(); ++ii)
- ret += XT::Common::VectorAbstraction<Vec>::get_entry(u, ii);
- return ret;
- }
-
- std::string short_id() const override final
- {
- return "hf";
- }
-
- std::string mn_name() const override final
- {
- return "hfm" + XT::Common::to_string(dimRange);
- }
-
- std::string pn_name() const override final
- {
- return "hfp" + XT::Common::to_string(dimRange);
- }
-
// get indices of all faces that contain point v
std::vector<size_t> get_face_indices(const DomainType& v) const
{
@@ -537,33 +572,8 @@ public:
using BaseType::u_iso;
- template <class Vec>
- std::enable_if_t<XT::Common::is_vector<Vec>::value, void> u_iso(Vec& ret) const
- {
- auto ret_range = u_iso();
- using V = XT::Common::VectorAbstraction<Vec>;
- for (size_t ii = 0; ii < dimRange; ++ii)
- V::set_entry(ret, ii, ret_range[ii]);
- }
-
- void ensure_min_density(DynamicRangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t ii = 0; ii < dimRange; ++ii)
- if (u[ii] < u_iso_min[ii])
- u[ii] = u_iso_min[ii];
- }
-
- void ensure_min_density(RangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t ii = 0; ii < dimRange; ++ii)
- if (u[ii] < u_iso_min[ii])
- u[ii] = u_iso_min[ii];
- }
-
- virtual DynamicRangeType
- get_moment_vector(const std::function<RangeFieldType(DomainType, bool)>& psi) const override final
+ DynamicRangeType
+ get_u(const std::function<RangeFieldType(DomainType, std::array<int, dimDomain>)>& psi) const override final
{
DynamicRangeType ret(dimRange, 0.);
const auto merged_quads = XT::Data::merged_quadrature(quadratures_);
@@ -573,13 +583,18 @@ public:
const auto& quad_point = *it;
const auto& v = quad_point.position();
const auto val = evaluate_on_face(v, face_index);
- const auto factor = psi(v, is_negative(it)) * quad_point.weight();
+ const auto factor = psi(v, has_fixed_sign(it)) * quad_point.weight();
for (size_t ii = 0; ii < 3; ++ii)
ret[vertices[ii]->index()] += val[ii] * factor;
}
return ret;
}
+ std::array<int, dimDomain> has_fixed_sign(const size_t index) const override final
+ {
+ return BaseType::triangle_has_fixed_sign(index);
+ }
+
protected:
template <class VertexVectorType>
bool calculate_barycentric_coordinates(const DomainType& v, const VertexVectorType& vertices, DomainType& ret) const
@@ -653,10 +668,10 @@ protected:
return ret;
}
- virtual void parallel_quadrature(const QuadraturesType& quadratures,
- MatrixType& matrix,
- const size_t v_index,
- const bool reflecting = false) const override final
+ void parallel_quadrature(const QuadraturesType& quadratures,
+ MatrixType& matrix,
+ const size_t v_index,
+ const bool reflecting = false) const override final
{
const auto& faces = triangulation_.faces();
size_t num_threads = std::min(XT::Common::threadManager().max_threads(), faces.size());
@@ -689,12 +704,12 @@ protected:
} // threads
} // void parallel_quadrature(...)
- virtual void calculate_in_thread_hat(std::vector<LocalMatrixType>& local_matrices,
- const QuadraturesType& quadratures,
- const size_t v_index,
- const std::vector<size_t>& decomposition,
- const size_t ii,
- const bool reflecting) const
+ void calculate_in_thread_hat(std::vector<LocalMatrixType>& local_matrices,
+ const QuadraturesType& quadratures,
+ const size_t v_index,
+ const std::vector<size_t>& decomposition,
+ const size_t ii,
+ const bool reflecting) const
{
const auto& reflected_indices = triangulation_.reflected_face_indices();
for (size_t face_index = decomposition[ii]; face_index < decomposition[ii + 1]; ++face_index) {
@@ -718,9 +733,9 @@ protected:
} // faces
} // void calculate_in_thread(...)
- virtual void integrated_initializer_thread(DynamicRangeType& local_range,
- const std::vector<MergedQuadratureIterator>& decomposition,
- const size_t ii) const override final
+ void integrated_initializer_thread(DynamicRangeType& local_range,
+ const std::vector<MergedQuadratureIterator>& decomposition,
+ const size_t ii) const override final
{
const auto& faces = triangulation_.faces();
for (auto it = decomposition[ii]; it != decomposition[ii + 1]; ++it) {
diff --git a/dune/gdt/test/momentmodels/basisfunctions/interface.hh b/dune/gdt/test/momentmodels/basisfunctions/interface.hh
index 61ee3c8b1d2b8fad9ce751d2c849f712aac69a5d..587de6be185f332b721d423187cf227297665334 100644
--- a/dune/gdt/test/momentmodels/basisfunctions/interface.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/interface.hh
@@ -214,7 +214,7 @@ public:
}
// returns <b>, where b is the basis functions vector
- virtual DynamicRangeType integrated() const
+ virtual const DynamicRangeType& integrated() const
{
return integrated_;
}
@@ -233,21 +233,30 @@ public:
return ret;
}
- virtual DynamicRangeType get_moment_vector(const std::function<RangeFieldType(DomainType, bool)>& psi) const
+ // Get moment vector from distribution function psi. psi also takes an array of bools as argument to decide which
+ // value to use at the discontinuities. Is currently ignored in all testcases except the Heaviside test in one
+ // dimension, so we only have to test if we are on the positive or negative interval touching 0. For more general
+ // testcases, this would have to be adapted.
+ virtual DynamicRangeType get_u(const std::function<RangeFieldType(DomainType, std::array<int, dimDomain>)>& psi) const
{
DynamicRangeType ret(dimRange, 0.);
const auto merged_quads = XT::Data::merged_quadrature(quadratures());
for (auto it = merged_quads.begin(); it != merged_quads.end(); ++it) {
const auto& quad_point = *it;
const auto& v = quad_point.position();
- ret += evaluate(v, it.first_index()) * psi(v, is_negative(it)) * quad_point.weight();
+ ret += evaluate(v, it.first_index()) * psi(v, has_fixed_sign(it.first_index())) * quad_point.weight();
}
return ret;
}
- virtual bool is_negative(const MergedQuadratureIterator& /*it*/) const
+ // Tests (per coordinate direction) whether the positions in the quadrature with index index are all negative or
+ // positive. Returns -1 if the entries are all non-positive, 1 if the entries are all non-negative, and 0 else.
+ virtual std::array<int, dimDomain> has_fixed_sign(const size_t /*index*/) const
{
- return false;
+ std::array<int, dimDomain> ret;
+ for (size_t dd = 0; dd < dimDomain; ++dd)
+ ret[dd] = 0;
+ return ret;
}
virtual FieldVector<MatrixType, dimFlux> flux_matrix() const
@@ -347,6 +356,11 @@ public:
}
}
+ virtual bool adjust_alpha_to_ensure_min_density(RangeType& /*alpha*/, const RangeFieldType /*min_density*/) const
+ {
+ return false;
+ }
+
// Volume of integration domain. For the Mn models it is important that u_iso has density 1. If the basis is exactly
// integrated, we thus use the exact unit ball volume. If the basis is only integrated by quadrature, we have to use
// <1> as volume to get a density of 1.
@@ -355,7 +369,7 @@ public:
return unit_ball_volume_exact();
}
- virtual DynamicRangeType u_iso() const
+ virtual const DynamicRangeType& u_iso() const
{
return u_iso_;
}
diff --git a/dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh b/dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh
index 2c5b7ff8c530e65ad6532ce8c8dc74dc7ab37e8b..20a69bc3a8a6769ee9a12b2d5d9ba55d1d53ee23 100644
--- a/dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh
@@ -13,6 +13,7 @@
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/null.hpp>
+#include <boost/math/special_functions/lambert_w.hpp>
#if HAVE_QHULL
# include <dune/xt/common/disable_warnings.hh>
@@ -29,6 +30,200 @@ namespace Dune {
namespace GDT {
+template <class DomainFieldType,
+ size_t domainDim,
+ class RangeFieldType,
+ size_t rangeDim,
+ size_t dimFlux,
+ EntropyType entropy,
+ size_t sizeBlock>
+class PartialMomentBasisBase
+ : public MomentBasisInterface<DomainFieldType, domainDim, RangeFieldType, rangeDim, 1, dimFlux, entropy>
+{
+ using BaseType = MomentBasisInterface<DomainFieldType, domainDim, RangeFieldType, rangeDim, 1, dimFlux, entropy>;
+ using ThisType = PartialMomentBasisBase;
+
+public:
+ using BaseType::dimDomain;
+ using BaseType::dimRange;
+ static constexpr size_t block_size = sizeBlock;
+ static constexpr size_t num_blocks = dimRange / block_size;
+
+ using typename BaseType::DomainType;
+ using typename BaseType::DynamicRangeType;
+ using typename BaseType::RangeType;
+ using typename BaseType::StringifierType;
+ using LocalVectorType = XT::Common::FieldVector<RangeFieldType, block_size>;
+
+ template <class... Args>
+ PartialMomentBasisBase(Args&&... args)
+ : BaseType(std::forward<Args>(args)...)
+ {}
+
+ // evaluate on spherical triangle face_index
+ DynamicRangeType evaluate(const DomainType& v, const size_t face_index) const override final
+ {
+ DynamicRangeType ret(dimRange, 0);
+ const auto local_eval = evaluate_on_face(v, face_index);
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[block_size * face_index + ii] = local_eval[ii];
+ return ret;
+ } // ... evaluate(...)
+
+ // evaluate on spherical triangle face_index
+ LocalVectorType evaluate_on_face(const DomainType& v, const size_t /*face_index*/) const
+ {
+ LocalVectorType ret;
+ ret[0] = 1.;
+ for (size_t ii = 1; ii < block_size; ++ii)
+ ret[ii] = v[ii - 1];
+ return ret;
+ } // ... evaluate(...)
+
+ static StringifierType stringifier()
+ {
+ return [](const DynamicRangeType& val) {
+ RangeFieldType psi(0);
+ for (size_t ii = 0; ii < dimRange; ii += block_size)
+ psi += val[ii];
+ return XT::Common::to_string(psi, 15);
+ };
+ } // ... stringifier()
+
+ DynamicRangeType alpha_one() const override final
+ {
+ DynamicRangeType ret(dimRange, 0.);
+ for (size_t ii = 0; ii < dimRange; ii += block_size)
+ ret[ii] = 1.;
+ return ret;
+ }
+
+ RangeFieldType density(const RangeType& u) const override final
+ {
+ RangeFieldType ret(0.);
+ for (size_t ii = 0; ii < dimRange; ii += block_size)
+ ret += u[ii];
+ return ret;
+ }
+
+ RangeFieldType density(const DynamicRangeType& u) const override final
+ {
+ RangeFieldType ret(0.);
+ for (size_t ii = 0; ii < dimRange; ii += block_size)
+ ret += u[ii];
+ return ret;
+ }
+
+ RangeFieldType density(const XT::Common::BlockedFieldVector<RangeFieldType, num_blocks, block_size>& u) const
+ {
+ RangeFieldType ret(0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ ret += u.block(jj)[0];
+ return ret;
+ }
+
+ RangeFieldType min_density(const XT::Common::BlockedFieldVector<RangeFieldType, num_blocks, block_size>& u) const
+ {
+ RangeFieldType ret(u.block(0)[0]);
+ for (size_t jj = 1; jj < num_blocks; ++jj)
+ ret = std::min(ret, u.block(jj)[0]);
+ return ret;
+ }
+
+ using BaseType::u_iso;
+
+ // For the partial moments, we might not be able to solve the optimization problem for some moments where the density
+ // on one interval/spherical triangle is very low. The overall density might be much higher than the density on that
+ // triangle, so we specialize this function.
+ void ensure_min_density(DynamicRangeType& u, const RangeFieldType min_density) const override final
+ {
+ const auto u_iso_min = u_iso() * min_density;
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto block_density = u[block_size * jj];
+ if (block_density < u_iso_min[block_size * jj]) {
+ for (size_t ii = 0; ii < block_size; ++ii)
+ u[block_size * jj + ii] = u_iso_min[block_size * jj + ii];
+ }
+ }
+ }
+
+ // For the partial moments, we might not be able to solve the optimization problem for some moments where the density
+ // on one interval/spherical triangle is very low. The overall density might be much higher than the density on that
+ // triangle, so we specialize this function.
+ void ensure_min_density(RangeType& u, const RangeFieldType min_density) const override final
+ {
+ const auto u_iso_min = u_iso() * min_density;
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto block_density = u[block_size * jj];
+ if (block_density < u_iso_min[block_size * jj]) {
+ for (size_t ii = 0; ii < block_size; ++ii)
+ u[block_size * jj + ii] = u_iso_min[block_size * jj + ii];
+ }
+ }
+ }
+
+ using BaseType::has_fixed_sign;
+
+ DynamicRangeType
+ get_u(const std::function<RangeFieldType(DomainType, std::array<int, dimDomain>)>& psi) const override final
+ {
+ DynamicRangeType ret(dimRange, 0);
+ for (size_t jj = 0; jj < quadratures_.size(); ++jj) {
+ const size_t offset = jj * block_size;
+ for (const auto& quad_point : quadratures_[jj]) {
+ const auto& v = quad_point.position();
+ const auto basis_val = evaluate_on_face(v, jj);
+ const auto psi_val = psi(v, has_fixed_sign(jj));
+ const auto factor = psi_val * quad_point.weight();
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[offset + ii] += basis_val[ii] * factor;
+ }
+ }
+ return ret;
+ }
+
+ std::string short_id() const override final
+ {
+ return "pm";
+ }
+
+ std::string mn_name() const override final
+ {
+ return "pmm" + XT::Common::to_string(dimRange);
+ }
+
+ std::string pn_name() const override final
+ {
+ return "pmp" + XT::Common::to_string(dimRange);
+ }
+
+protected:
+ using BaseType::quadratures_;
+}; // class PartialMomentBasisBase<...>
+
+template <class DomainFieldType,
+ size_t domainDim,
+ class RangeFieldType,
+ size_t rangeDim,
+ size_t dimFlux,
+ EntropyType entropy,
+ size_t sizeBlock>
+constexpr size_t
+ PartialMomentBasisBase<DomainFieldType, domainDim, RangeFieldType, rangeDim, dimFlux, entropy, sizeBlock>::
+ num_blocks;
+
+template <class DomainFieldType,
+ size_t domainDim,
+ class RangeFieldType,
+ size_t rangeDim,
+ size_t dimFlux,
+ EntropyType entropy,
+ size_t sizeBlock>
+constexpr size_t
+ PartialMomentBasisBase<DomainFieldType, domainDim, RangeFieldType, rangeDim, dimFlux, entropy, sizeBlock>::
+ block_size;
+
+
template <class DomainFieldType,
size_t dimDomain,
class RangeFieldType,
@@ -49,20 +244,16 @@ template <class DomainFieldType,
size_t dimFlux,
EntropyType entropy>
class PartialMomentBasis<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, dimFlux, 1, entropy>
- : public MomentBasisInterface<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, dimFlux, entropy>
+ : public PartialMomentBasisBase<DomainFieldType, 1, RangeFieldType, rangeDim, dimFlux, entropy, 2>
{
+ using BaseType = PartialMomentBasisBase<DomainFieldType, 1, RangeFieldType, rangeDim, dimFlux, entropy, 2>;
+
public:
- static constexpr size_t dimDomain = 1;
- static constexpr size_t dimRange = rangeDim;
- static constexpr size_t dimRangeCols = rangeDimCols;
+ using BaseType::dimDomain;
+ using BaseType::dimRange;
static_assert(!(dimRange % 2), "dimRange has to be even!");
- static constexpr size_t num_intervals = dimRange / 2;
-
-private:
- using BaseType =
- MomentBasisInterface<DomainFieldType, dimDomain, RangeFieldType, dimRange, dimRangeCols, dimFlux, entropy>;
+ static constexpr size_t num_intervals = BaseType::num_blocks;
-public:
using typename BaseType::DomainType;
using typename BaseType::DynamicRangeType;
using typename BaseType::MatrixType;
@@ -70,15 +261,9 @@ public:
using typename BaseType::RangeType;
using typename BaseType::SphericalTriangulationType;
using typename BaseType::StringifierType;
- using typename BaseType::VisualizerType;
using LocalVectorType = FieldVector<RangeFieldType, 2>;
using PartitioningType = typename BaseType::Partitioning1dType;
- static std::string static_id()
- {
- return "pcw";
- }
-
static size_t default_quad_order()
{
return 15;
@@ -105,6 +290,8 @@ public:
BaseType::initialize_base_values();
}
+ using BaseType::evaluate;
+
DynamicRangeType evaluate(const DomainType& v) const override final
{
for (size_t ii = 0; ii < num_intervals; ++ii)
@@ -114,41 +301,54 @@ public:
return DynamicRangeType();
} // ... evaluate(...)
- // evaluate on interval ii
- DynamicRangeType evaluate(const DomainType& v, const size_t ii) const override final
- {
- DynamicRangeType ret(dimRange, 0);
- const auto local_ret = evaluate_on_interval(v, ii);
- ret[2 * ii] = local_ret[0];
- ret[2 * ii + 1] = local_ret[1];
- return ret;
- } // ... evaluate(...)
-
- LocalVectorType evaluate_on_interval(const DomainType& v, const size_t /*ii*/) const
+ bool adjust_alpha_to_ensure_min_density(RangeType& alpha, const RangeFieldType psi_min) const override final
{
- return LocalVectorType{{1, v[0]}};
- } // ... evaluate(...)
-
- LocalVectorType evaluate_on_face(const DomainType& v, const size_t ii) const
- {
- return evaluate_on_interval(v, ii);
- } // ... evaluate(...)
-
- DynamicRangeType integrated() const override final
- {
- DynamicRangeType ret(dimRange, 0);
+ bool changed = false;
+ const auto alpha_min = std::log(psi_min);
for (size_t ii = 0; ii < num_intervals; ++ii) {
- ret[2 * ii] = partitioning_[ii + 1] - partitioning_[ii];
- ret[2 * ii + 1] = (std::pow(partitioning_[ii + 1], 2) - std::pow(partitioning_[ii], 2)) / 2.;
- }
- return ret;
- }
-
- virtual bool
- is_negative(const typename XT::Data::MergedQuadrature<RangeFieldType, dimDomain>::MergedQuadratureIterator& it)
- const override final
- {
- return it.first_index() < num_intervals / 2;
+ auto& alpha_0 = alpha[2 * ii];
+ auto& alpha_1 = alpha[2 * ii + 1];
+ const auto mu_i = partitioning_[ii];
+ const auto mu_ip1 = partitioning_[ii + 1];
+ const auto alpha_left = alpha_0 + mu_i * alpha_1;
+ const auto alpha_right = alpha_0 + mu_ip1 * alpha_1;
+ const bool min_is_left = alpha_left < alpha_right;
+ const auto alpha_min_ii = min_is_left ? alpha_left : alpha_right;
+ const auto alpha_max_ii = min_is_left ? alpha_right : alpha_left;
+ if (alpha_min_ii < alpha_min) {
+ if (XT::Common::FloatCmp::le(alpha_max_ii, alpha_min)) {
+ alpha_0 = alpha_min;
+ alpha_1 = 0.;
+ changed = true;
+ } else {
+ // We know that alpha_1 != 0 because alpha_max_ii > alpha_min and alpha_min_ii < alpha_min
+ const auto rho_ii = -(std::exp(alpha_0 + mu_i * alpha_1) - std::exp(alpha_0 + mu_ip1 * alpha_1)) / alpha_1;
+ if (std::isnan(rho_ii) || std::isinf(rho_ii))
+ DUNE_THROW(Dune::MathError, "Inf or nan in rho!");
+ const auto h = (mu_ip1 - mu_i);
+ if (XT::Common::FloatCmp::le(rho_ii, psi_min * h)) {
+ alpha_0 = alpha_min;
+ alpha_1 = 0.;
+ changed = true;
+ continue;
+ }
+ // get positive slope
+#if 0
+ // Set minimum to alpha_min, leave max alpha unchanged
+ alpha_1 = (alpha_max_ii - alpha_min_ii) / h;
+#else
+ // Set minimum to alpha_min, leave rho unchanged
+ alpha_1 = -1. / h * boost::math::lambert_wm1(-h * psi_min / rho_ii * std::exp(-h * psi_min / rho_ii))
+ - psi_min / rho_ii;
+#endif
+ if (!min_is_left)
+ alpha_1 *= -1.; // slope has to be negative in this case
+ alpha_0 = alpha_min - (min_is_left ? mu_i : mu_ip1) * alpha_1;
+ changed = true;
+ }
+ }
+ } // ii
+ return changed;
}
// returns matrix with entries <h_i h_j>
@@ -250,105 +450,11 @@ public:
return ret;
}
- static StringifierType stringifier()
- {
- return [](const RangeType& val) {
- RangeFieldType psi(0);
- for (size_t ii = 0; ii < dimRange; ii += 2)
- psi += val[ii];
- return XT::Common::to_string(psi, 15);
- };
- } // ... stringifier()
-
const PartitioningType& partitioning() const
{
return partitioning_;
}
- DynamicRangeType alpha_one() const override final
- {
- DynamicRangeType ret(dimRange, 0);
- for (size_t ii = 0; ii < dimRange; ii += 2)
- ret[ii] = 1.;
- return ret;
- }
-
- RangeFieldType density(const DynamicRangeType& u) const override final
- {
- RangeFieldType ret(0.);
- for (size_t ii = 0; ii < dimRange; ii += 2) {
- ret += u[ii];
- }
- return ret;
- }
-
- RangeFieldType density(const RangeType& u) const
- {
- RangeFieldType ret(0.);
- for (size_t ii = 0; ii < dimRange; ii += 2) {
- ret += u[ii];
- }
- return ret;
- }
-
- RangeFieldType density(const XT::Common::BlockedFieldVector<RangeFieldType, num_intervals, 2>& u) const
- {
- RangeFieldType ret(0.);
- for (size_t jj = 0; jj < num_intervals; ++jj)
- ret += u.block(jj)[0];
- return ret;
- }
-
- RangeFieldType min_density(const XT::Common::BlockedFieldVector<RangeFieldType, num_intervals, 2>& u) const
- {
- RangeFieldType ret(u.block(0)[0]);
- for (size_t jj = 1; jj < num_intervals; ++jj)
- ret = std::min(ret, u.block(jj)[0]);
- return ret;
- }
-
- using BaseType::u_iso;
-
- // For the partial moments, we might not be able to solve the optimization problem for some moments where the density
- // on one interval/spherical triangle is very low. The overall density might be much higher than the density on that
- // triangle, so we specialize this function.
- void ensure_min_density(DynamicRangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t jj = 0; jj < num_intervals; ++jj) {
- if (u[2 * jj] < u_iso_min[2 * jj]) {
- u[2 * jj] = u_iso_min[2 * jj];
- u[2 * jj + 1] = u_iso_min[2 * jj + 1];
- }
- }
- }
-
- void ensure_min_density(RangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t jj = 0; jj < num_intervals; ++jj) {
- if (u[2 * jj] < u_iso_min[2 * jj]) {
- u[2 * jj] = u_iso_min[2 * jj];
- u[2 * jj + 1] = u_iso_min[2 * jj + 1];
- }
- }
- }
-
- std::string short_id() const override final
- {
- return "pm";
- }
-
- std::string mn_name() const override final
- {
- return "pmm" + XT::Common::to_string(dimRange);
- }
-
- std::string pn_name() const override final
- {
- return "pmp" + XT::Common::to_string(dimRange);
- }
-
// get indices of all faces that contain point
std::vector<size_t> get_face_indices(const DomainType& v) const
{
@@ -360,6 +466,21 @@ public:
return face_indices;
}
+ DynamicRangeType integrated_initializer(const QuadraturesType& /*quadratures*/) const override final
+ {
+ DynamicRangeType ret(dimRange, 0);
+ for (size_t ii = 0; ii < num_intervals; ++ii) {
+ ret[2 * ii] = partitioning_[ii + 1] - partitioning_[ii];
+ ret[2 * ii + 1] = (std::pow(partitioning_[ii + 1], 2) - std::pow(partitioning_[ii], 2)) / 2.;
+ }
+ return ret;
+ }
+
+ std::array<int, dimDomain> has_fixed_sign(const size_t index) const override final
+ {
+ return BaseType::interval_has_fixed_sign(index, num_intervals);
+ }
+
private:
const PartitioningType partitioning_;
using BaseType::quadratures_;
@@ -372,32 +493,34 @@ template <class DomainFieldType,
size_t dimFlux,
EntropyType entropy>
constexpr size_t
- PartialMomentBasis<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, dimFlux, 1, entropy>::dimRange;
+ PartialMomentBasis<DomainFieldType, 1, RangeFieldType, rangeDim, rangeDimCols, dimFlux, 1, entropy>::num_intervals;
template <class DomainFieldType, class RangeFieldType, size_t refinements, size_t dimFlux, EntropyType entropy>
class PartialMomentBasis<DomainFieldType, 3, RangeFieldType, refinements, 1, dimFlux, 1, entropy>
- : public MomentBasisInterface<DomainFieldType,
- 3,
- RangeFieldType,
- OctaederStatistics<refinements>::num_faces() * 4,
- 1,
- dimFlux,
- entropy>
+ : public PartialMomentBasisBase<DomainFieldType,
+ 3,
+ RangeFieldType,
+ OctaederStatistics<refinements>::num_faces() * 4,
+ dimFlux,
+ entropy,
+ 4>
{
-public:
- static const size_t dimDomain = 3;
- static const size_t dimRange = OctaederStatistics<refinements>::num_faces() * 4;
- static const size_t dimRangeCols = 1;
- static constexpr size_t block_size = 4;
- static constexpr size_t num_blocks = dimRange / block_size;
- static constexpr size_t num_refinements = refinements;
-
-private:
- using BaseType =
- MomentBasisInterface<DomainFieldType, dimDomain, RangeFieldType, dimRange, dimRangeCols, dimFlux, entropy>;
+ using BaseType = PartialMomentBasisBase<DomainFieldType,
+ 3,
+ RangeFieldType,
+ OctaederStatistics<refinements>::num_faces() * 4,
+ dimFlux,
+ entropy,
+ 4>;
using ThisType = PartialMomentBasis;
public:
+ using BaseType::block_size;
+ using BaseType::dimDomain;
+ using BaseType::dimRange;
+ using BaseType::num_blocks;
+ static constexpr size_t num_refinements = refinements;
+
using TriangulationType = typename BaseType::SphericalTriangulationType;
using typename BaseType::DomainType;
using typename BaseType::DynamicRangeType;
@@ -406,7 +529,6 @@ public:
using typename BaseType::QuadratureType;
using typename BaseType::RangeType;
using typename BaseType::StringifierType;
- using typename BaseType::VisualizerType;
using BlockRangeType = XT::Common::FieldVector<RangeFieldType, block_size>;
using BlockPlaneCoefficientsType = typename std::vector<std::pair<BlockRangeType, RangeFieldType>>;
using PlaneCoefficientsType = XT::Common::FieldVector<BlockPlaneCoefficientsType, num_blocks>;
@@ -456,6 +578,8 @@ public:
BaseType::initialize_base_values();
}
+ using BaseType::evaluate;
+
DynamicRangeType evaluate(const DomainType& v) const override final
{
DynamicRangeType ret(dimRange, 0);
@@ -467,129 +591,11 @@ public:
return evaluate(v, face_indices[0]);
} // ... evaluate(...)
- // evaluate on spherical triangle face_index
- DynamicRangeType evaluate(const DomainType& v, const size_t face_index) const override final
- {
- DynamicRangeType ret(dimRange, 0);
- const auto local_eval = evaluate_on_face(v, face_index);
- assert(4 * face_index + 3 < ret.size());
- for (size_t ii = 0; ii < 4; ++ii)
- ret[4 * face_index + ii] = local_eval[ii];
- return ret;
- } // ... evaluate(...)
-
- // evaluate on spherical triangle face_index
- LocalVectorType evaluate_on_face(const DomainType& v, const size_t /*face_index*/) const
- {
- LocalVectorType ret;
- ret[0] = 1.;
- for (size_t ii = 1; ii < 4; ++ii)
- ret[ii] = v[ii - 1];
- return ret;
- } // ... evaluate(...)
-
- static StringifierType stringifier()
- {
- return [](const DynamicRangeType& val) {
- RangeFieldType psi(0);
- for (size_t ii = 0; ii < dimRange; ii += 4)
- psi += val[ii];
- return XT::Common::to_string(psi, 15);
- };
- } // ... stringifier()
-
RangeFieldType unit_ball_volume() const override final
{
return BaseType::unit_ball_volume_quad();
}
- DynamicRangeType alpha_one() const override final
- {
- DynamicRangeType ret(dimRange, 0.);
- for (size_t ii = 0; ii < dimRange; ii += 4)
- ret[ii] = 1.;
- return ret;
- }
-
- virtual RangeFieldType density(const RangeType& u) const
- {
- RangeFieldType ret(0.);
- for (size_t ii = 0; ii < dimRange; ii += 4)
- ret += u[ii];
- return ret;
- }
-
- RangeFieldType density(const DynamicRangeType& u) const override final
- {
- RangeFieldType ret(0.);
- for (size_t ii = 0; ii < dimRange; ii += 4)
- ret += u[ii];
- return ret;
- }
-
- RangeFieldType density(const XT::Common::BlockedFieldVector<RangeFieldType, dimRange / 4, 4>& u) const
- {
- RangeFieldType ret(0.);
- for (size_t jj = 0; jj < dimRange / 4; ++jj)
- ret += u.block(jj)[0];
- return ret;
- }
-
- RangeFieldType min_density(const XT::Common::BlockedFieldVector<RangeFieldType, dimRange / 4, 4>& u) const
- {
- RangeFieldType ret(u.block(0)[0]);
- for (size_t jj = 1; jj < dimRange / 4; ++jj)
- ret = std::min(ret, u.block(jj)[0]);
- return ret;
- }
-
- using BaseType::u_iso;
-
- // For the partial moments, we might not be able to solve the optimization problem for some moments where the density
- // on one interval/spherical triangle is very low. The overall density might be much higher than the density on that
- // triangle, so we specialize this function.
- void ensure_min_density(DynamicRangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- const auto block_density = u[4 * jj];
- if (block_density < u_iso_min[4 * jj]) {
- for (size_t ii = 0; ii < block_size; ++ii)
- u[4 * jj + ii] = u_iso_min[4 * jj + ii];
- }
- }
- }
-
- // For the partial moments, we might not be able to solve the optimization problem for some moments where the density
- // on one interval/spherical triangle is very low. The overall density might be much higher than the density on that
- // triangle, so we specialize this function.
- void ensure_min_density(RangeType& u, const RangeFieldType min_density) const override final
- {
- const auto u_iso_min = u_iso() * min_density;
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- const auto block_density = u[4 * jj];
- if (block_density < u_iso_min[4 * jj]) {
- for (size_t ii = 0; ii < block_size; ++ii)
- u[4 * jj + ii] = u_iso_min[4 * jj + ii];
- }
- }
- }
-
- std::string short_id() const override final
- {
- return "pm";
- }
-
- std::string mn_name() const override final
- {
- return "pmm" + XT::Common::to_string(dimRange);
- }
-
- std::string pn_name() const override final
- {
- return "pmp" + XT::Common::to_string(dimRange);
- }
-
std::vector<size_t> get_face_indices(const DomainType& v) const
{
return triangulation_.get_face_indices(v);
@@ -637,24 +643,6 @@ public:
threads[jj].join();
}
- virtual DynamicRangeType
- get_moment_vector(const std::function<RangeFieldType(DomainType, bool)>& psi) const override final
- {
- DynamicRangeType ret(dimRange, 0);
- for (size_t jj = 0; jj < quadratures_.size(); ++jj) {
- const size_t offset = jj * block_size;
- for (const auto& quad_point : quadratures_[jj]) {
- const auto& v = quad_point.position();
- const auto basis_val = evaluate_on_face(v, jj);
- const auto psi_val = psi(v, false);
- const auto factor = psi_val * quad_point.weight();
- for (size_t ii = 0; ii < block_size; ++ii)
- ret[offset + ii] += basis_val[ii] * factor;
- }
- }
- return ret;
- }
-
std::unique_ptr<BlockMatrixType> block_mass_matrix() const
{
auto block_matrix = std::make_unique<BlockMatrixType>();
@@ -697,6 +685,10 @@ public:
return B_kinetic;
} // ... kinetic_flux_matrices()
+ std::array<int, dimDomain> has_fixed_sign(const size_t index) const override final
+ {
+ return BaseType::triangle_has_fixed_sign(index);
+ }
private:
void calculate_plane_coefficients_block(std::vector<XT::Common::FieldVector<RangeFieldType, block_size>>& points,
@@ -821,13 +813,6 @@ private:
mutable PlaneCoefficientsType plane_coefficients_;
}; // class PartialMomentBasis<DomainFieldType, 3, ...>
-template <class DomainFieldType, class RangeFieldType, size_t refinements, size_t dimFlux, EntropyType entropy>
-constexpr size_t PartialMomentBasis<DomainFieldType, 3, RangeFieldType, refinements, 1, dimFlux, 1, entropy>::dimRange;
-
-template <class DomainFieldType, class RangeFieldType, size_t refinements, size_t dimFlux, EntropyType entropy>
-constexpr size_t
- PartialMomentBasis<DomainFieldType, 3, RangeFieldType, refinements, 1, dimFlux, 1, entropy>::num_blocks;
-
template <class DomainFieldType, class RangeFieldType, size_t refinements, size_t dimFlux, EntropyType entropy>
constexpr size_t
PartialMomentBasis<DomainFieldType, 3, RangeFieldType, refinements, 1, dimFlux, 1, entropy>::num_refinements;
diff --git a/dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh b/dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh
index bdd0afa436fecd26403386ce6962e7bfb8b0a458..4a4f5e9cf5d0e4741a94864a8ac36c1e1bb5a19d 100644
--- a/dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh
@@ -105,13 +105,6 @@ public:
return ret;
} // ... evaluate(...)
- DynamicRangeType integrated() const override final
- {
- DynamicRangeType ret(dimRange, 0.);
- ret[0] = std::sqrt(4. * M_PI);
- return ret;
- }
-
MatrixType mass_matrix() const override
{
MatrixType M(dimRange, dimRange, 0);
@@ -167,6 +160,14 @@ public:
return "p" + XT::Common::to_string(order);
}
+
+ DynamicRangeType integrated_initializer() const override final
+ {
+ DynamicRangeType ret(dimRange, 0.);
+ ret[0] = std::sqrt(4. * M_PI);
+ return ret;
+ }
+
private:
static RangeFieldType A_lm(const int l, const int m)
{
diff --git a/dune/gdt/test/momentmodels/kinetictransport/base.hh b/dune/gdt/test/momentmodels/kinetictransport/base.hh
index c9705f94278d971d8877fdb332d2d2d35afe3073..5fa1e676b1ffb1e12946eb3b4de44fdc0a9aac8e 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/base.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/base.hh
@@ -41,6 +41,7 @@ public:
using typename BaseType::BasisDomainType;
using typename BaseType::DomainFieldType;
using typename BaseType::DomainType;
+ using typename BaseType::DynamicRangeType;
using typename BaseType::GenericFluxFunctionType;
using typename BaseType::GenericFunctionType;
using typename BaseType::MatrixType;
@@ -162,19 +163,19 @@ public:
// Thus, the initial value of the n-th moment is basis_integrated * psi_vac.
std::unique_ptr<InitialValueType> initial_values() const override
{
- RangeReturnType value = basis_functions_.integrated() * psi_vac_;
+ const auto value = basis_functions_.integrated() * psi_vac_;
return std::make_unique<GenericFunctionType>(
[](const XT::Common::Parameter&) { return 0; },
- [value](const DomainType&, const XT::Common::Parameter&) { return value; });
+ [value](const DomainType&, DynamicRangeType& ret, const XT::Common::Parameter&) { ret = value; });
} // ... initial_values()
// Use a constant vacuum concentration basis_integrated * psi_vac as default boundary value
std::unique_ptr<BoundaryValueType> boundary_values() const override
{
- RangeReturnType value = basis_functions_.integrated() * psi_vac_;
+ const auto value = basis_functions_.integrated() * psi_vac_;
return std::make_unique<GenericFunctionType>(
[](const XT::Common::Parameter&) { return 0; },
- [=](const DomainType&, const XT::Common::Parameter&) { return value; });
+ [=](const DomainType&, DynamicRangeType& ret, const XT::Common::Parameter&) { ret = value; });
} // ... boundary_values()
virtual BoundaryDistributionType boundary_distribution() const
diff --git a/dune/gdt/test/momentmodels/kinetictransport/planesource.hh b/dune/gdt/test/momentmodels/kinetictransport/planesource.hh
index 01906d8e43f53de2a98da736734bd35cc72f4377..715cbfbfbe977b0d13ea66f3091c43a6979bdcbf 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/planesource.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/planesource.hh
@@ -28,6 +28,7 @@ public:
using BaseType::dimRange;
using typename BaseType::ConstantScalarFunctionType;
using typename BaseType::DomainType;
+ using typename BaseType::DynamicRangeType;
using typename BaseType::GenericFunctionType;
using typename BaseType::InitialValueType;
using typename BaseType::MomentBasis;
@@ -66,18 +67,17 @@ public:
const size_t num_elements = XT::Common::from_string<std::vector<size_t>>(grid_cfg_["num_elements"])[0];
const RangeFieldType len_domain = upper_right[0] - lower_left[0];
const RangeFieldType vol_entity = len_domain / num_elements;
- RangeReturnType basis_integrated = basis_functions_.integrated();
+ const auto basis_integrated = basis_functions_.integrated();
const RangeFieldType domain_center = lower_left[0] + len_domain / 2;
// approximate delta function by constant value of 1/(2*vol_entity) on cells on both side of 0.
- const auto eval_func = [=](const DomainType& x, const XT::Common::Parameter&) {
- auto ret = basis_integrated;
+ const auto eval_func = [=](const DomainType& x, DynamicRangeType& ret, const XT::Common::Parameter&) {
+ ret = basis_integrated;
if (XT::Common::FloatCmp::ge(x[0], domain_center - vol_entity)
&& XT::Common::FloatCmp::le(x[0], domain_center + vol_entity))
ret *= psi_vac_ + 1. / (2. * vol_entity);
else
ret *= psi_vac_;
- return ret;
};
return std::make_unique<GenericFunctionType>(0, eval_func);
} // ... initial_values()
diff --git a/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh b/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
index 92e6da5972060d74592b08eba6187e622bdd3e23..a264c51b0125562666bf156c033a244cb0413500 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
@@ -82,16 +82,15 @@ public:
// at x = 3.
std::unique_ptr<BoundaryValueType> boundary_values() const override final
{
- return std::make_unique<GenericFunctionType>(1, [&](const DomainType& x, const XT::Common::Parameter&) {
- if (x[0] < 1.5) {
- static auto ret = helper<MomentBasis>::get_left_boundary_values(basis_functions_, psi_vac_, is_mn_model_);
- return ret;
- } else {
- auto ret = basis_functions_.integrated();
- ret *= psi_vac_;
- return ret;
- }
- });
+ return std::make_unique<GenericFunctionType>(
+ 1, [&](const DomainType& x, DynamicRangeType& ret, const XT::Common::Parameter&) {
+ if (x[0] < 1.5) {
+ helper<MomentBasis>::get_left_boundary_values(basis_functions_, psi_vac_, is_mn_model_, ret);
+ } else {
+ ret = basis_functions_.integrated();
+ ret *= psi_vac_;
+ }
+ });
} // ... boundary_values()
BoundaryDistributionType boundary_distribution() const override final
@@ -115,7 +114,9 @@ public:
RangeReturnType left_boundary_value() const
{
- return helper<MomentBasis>::get_left_boundary_values(basis_functions_, psi_vac_, is_mn_model_);
+ DynamicRangeType ret;
+ helper<MomentBasis>::get_left_boundary_values(basis_functions_, psi_vac_, is_mn_model_, ret);
+ return ret;
}
RangeFieldType t_end() const override
@@ -198,15 +199,16 @@ protected:
using helper_base::denominator;
using helper_base::g;
- static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
- const RangeFieldType& psi_vac,
- const bool is_mn_model)
+ static void get_left_boundary_values(const MomentBasisImp& basis_functions,
+ const RangeFieldType& psi_vac,
+ const bool is_mn_model,
+ DynamicRangeType& ret)
{
- DynamicRangeType ret(dimRange, 0);
// For the MN-Models, we have to use the quadrature also used in the optimization problem to guarantee
// realizability of the boundary_values.
// For the PN-Models, we do not have these issues and just use a very fine quadrature (which is not a performance
// problem as the integration is only done once).
+ std::fill(ret.begin(), ret.end(), 0.);
const auto& quadratures =
is_mn_model ? basis_functions.quadratures() : MomentBasisImp::gauss_lobatto_quadratures(100, 31);
for (size_t ii = 0; ii < quadratures.size(); ++ii) {
@@ -221,7 +223,6 @@ protected:
ret /= denominator();
// add small vacuum concentration to move away from realizable boundary
ret += basis_functions.integrated() * psi_vac;
- return ret;
}
static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& /*basis_functions*/,
@@ -246,11 +247,12 @@ protected:
using helper_base::integral_2;
using helper_base::integral_3;
- static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
- const RangeFieldType psi_vac,
- const bool /*is_mn_model*/)
+ static void get_left_boundary_values(const MomentBasisImp& basis_functions,
+ const RangeFieldType psi_vac,
+ const bool /*is_mn_model*/,
+ DynamicRangeType& ret)
{
- DynamicRangeType ret(dimRange, 0);
+ std::fill(ret.begin(), ret.end(), 0.);
for (size_t nn = 0; nn < dimRange; ++nn) {
const auto& partitioning = basis_functions.partitioning();
const auto vn = partitioning[nn];
@@ -265,7 +267,6 @@ protected:
}
// add small vacuum concentration to move away from realizable boundary
ret += basis_functions.integrated() * psi_vac;
- return ret;
} // ... get_left_boundary_values(...)
static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& basis_functions,
@@ -316,19 +317,19 @@ protected:
using helper_base::integral_2;
using helper_base::integral_3;
- static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
- const RangeFieldType psi_vac,
- const bool /*is_mn_model*/)
+ static void get_left_boundary_values(const MomentBasisImp& basis_functions,
+ const RangeFieldType psi_vac,
+ const bool /*is_mn_model*/,
+ DynamicRangeType& ret)
{
const auto& partitioning = basis_functions.partitioning();
- DynamicRangeType ret(dimRange, 0);
+ std::fill(ret.begin(), ret.end(), 0.);
for (size_t ii = 0; ii < dimRange / 2; ++ii) {
ret[2 * ii] = integral_1(partitioning[ii], partitioning[ii + 1]) / denominator();
ret[2 * ii + 1] = integral_2(partitioning[ii], partitioning[ii + 1]) / denominator();
}
// add small vacuum concentration to move away from realizable boundary
ret += basis_functions.integrated() * psi_vac;
- return ret;
}
static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& basis_functions,