diff --git a/dune/stuff/functions/spe10.hh b/dune/stuff/functions/spe10.hh
index 4b798b5fb0dca9b12e94c34037e9a32df2354afa..711bfdf0d7a9b6a0a5871d3c0e842b05cf92e39e 100644
--- a/dune/stuff/functions/spe10.hh
+++ b/dune/stuff/functions/spe10.hh
@@ -217,6 +217,84 @@ private:
}
}; // class Model1< ..., 2, ..., r, r >
+template <class EntityImp, class DomainFieldImp, size_t dim_domain, class RangeFieldImp, size_t r, size_t rC>
+class Model2 : public Stuff::GlobalFunctionInterface<EntityImp, DomainFieldImp, dim_domain, RangeFieldImp, r, rC>
+{
+ static_assert(r == rC, "");
+ static_assert(dim_domain == rC, "");
+ static_assert(dim_domain == 3, "");
+ typedef Stuff::GlobalFunctionInterface<EntityImp, DomainFieldImp, dim_domain, RangeFieldImp, r, rC> BaseType;
+
+public:
+ Model2(std::string data_filename = "perm_case2a.dat")
+ : deltas_{{6.096, 3.048, 0.6096}}
+ , permeability_(nullptr)
+ , permMatrix_(0.0)
+ , filename_(data_filename)
+ {
+ readPermeability();
+ }
+
+ virtual ~Model2()
+ {
+ delete permeability_;
+ permeability_ = nullptr;
+ }
+
+ //! currently used in gdt assembler
+ virtual void evaluate(const typename BaseType::DomainType& x,
+ typename BaseType::RangeType& diffusion) const final override
+ {
+
+ if (!permeability_) {
+ DSC_LOG_ERROR_0 << "The SPE10-permeability data file could not be opened. This file does\n"
+ << "not come with the dune-multiscale repository due to file size. To download it\n"
+ << "execute\n"
+ << "wget http://www.spe.org/web/csp/datasets/por_perm_case2a.zip\n"
+ << "unzip the file and move the file 'spe_perm.dat' to\n"
+ << "dune-multiscale/dune/multiscale/problems/spe10_permeability.dat!\n";
+ DUNE_THROW(IOError, "Data file for Groundwaterflow permeability could not be opened!");
+ }
+ // 3 is the maximum space dimension
+ for (size_t dim = 0; dim < dim_domain; ++dim)
+ permIntervalls_[dim] = std::floor(x[dim] / deltas_[dim]);
+
+ const int offset = permIntervalls_[0] + permIntervalls_[1] * 60 + permIntervalls_[2] * 220 * 60;
+ for (size_t dim = 0; dim < dim_domain; ++dim)
+ diffusion[dim][dim] = permeability_[offset + dim * 1122000];
+ }
+
+ virtual size_t order() const
+ {
+ return 0u;
+ }
+
+private:
+ void readPermeability()
+ {
+ std::ifstream file(filename_);
+ double val;
+ if (!file) { // file couldn't be opened
+ return;
+ }
+ file >> val;
+ int counter = 0;
+ permeability_ = new double[3366000];
+ while (!file.eof()) {
+ // keep reading until end-of-file
+ permeability_[counter++] = val;
+ file >> val; // sets EOF flag if no value found
+ }
+ file.close();
+ }
+
+ std::array<double, dim_domain> deltas_;
+ double* permeability_; //! TODO automatic memory
+ mutable typename BaseType::DomainType permIntervalls_;
+ mutable Dune::FieldMatrix<double, BaseType::DomainType::dimension, BaseType::DomainType::dimension> permMatrix_;
+ const std::string filename_;
+};
+
} // namespace Spe10
} // namespace Functions
diff --git a/dune/stuff/test/functions_spe10.cc b/dune/stuff/test/functions_spe10.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7be7c2c22fdcbedfee80b5199a02f32241e87d1e
--- /dev/null
+++ b/dune/stuff/test/functions_spe10.cc
@@ -0,0 +1,162 @@
+// This file is part of the dune-stuff project:
+// https://github.com/wwu-numerik/dune-stuff
+// Copyright holders: Rene Milk, Felix Schindler
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+#include "main.hxx"
+
+#include <memory>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/stuff/functions/interfaces.hh>
+#include <dune/stuff/functions/spe10.hh>
+
+// we need this nasty code generation because the testing::Types< ... > only accepts 50 arguments
+// and all combinations of functions and entities and dimensions and fieldtypes would be way too much
+#define TEST_STRUCT_GENERATOR(ftype, etype) \
+ template <class Spe10Model2FunctionType> \
+ struct ftype##etype##Test : public ::testing::Test \
+ { \
+ typedef typename Spe10Model2FunctionType::EntityType EntityType; \
+ typedef typename Spe10Model2FunctionType::LocalfunctionType LocalfunctionType; \
+ typedef typename Spe10Model2FunctionType::DomainFieldType DomainFieldType; \
+ static const size_t dimDomain = Spe10Model2FunctionType::dimDomain; \
+ typedef typename Spe10Model2FunctionType::DomainType DomainType; \
+ typedef typename Spe10Model2FunctionType::RangeFieldType RangeFieldType; \
+ static const size_t dimRange = Spe10Model2FunctionType::dimRange; \
+ static const size_t dimRangeCols = Spe10Model2FunctionType::dimRangeCols; \
+ typedef typename Spe10Model2FunctionType::RangeType RangeType; \
+ typedef typename Spe10Model2FunctionType::JacobianRangeType JacobianRangeType; \
+ \
+ void check() const \
+ { \
+ const Spe10Model2FunctionType zero; \
+ } \
+ };
+// TEST_STRUCT_GENERATOR
+
+
+#if HAVE_DUNE_GRID
+
+//# include <dune/grid/sgrid.hh>
+
+// typedef Dune::SGrid< 1, 1 >::Codim< 0 >::Entity DuneSGrid1dEntityType;
+// typedef Dune::SGrid< 2, 2 >::Codim< 0 >::Entity DuneSGrid2dEntityType;
+// typedef Dune::SGrid< 3, 3 >::Codim< 0 >::Entity DuneSGrid3dEntityType;
+
+// typedef testing::Types< Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 1, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 1, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 1, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 2, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 2, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 2, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 3, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 3, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid1dEntityType, double, 1, double, 3, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 1, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 1, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 1, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 2, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 2, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 2, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 3, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 3, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid2dEntityType, double, 2, double, 3, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 1, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 1, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 1, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 2, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 2, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 2, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 3, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 3, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneSGrid3dEntityType, double, 3, double, 3, 3 >
+// > Spe10Model2FunctionSGridEntityTypes;
+
+// TEST_STRUCT_GENERATOR(Spe10Model2FunctionTest, SGridEntity)
+// TYPED_TEST_CASE(Spe10Model2FunctionTestSGridEntityTest, Spe10Model2FunctionSGridEntityTypes);
+// TYPED_TEST(Spe10Model2FunctionTestSGridEntityTest, provides_required_methods) {
+// this->check();
+//}
+
+#include <dune/grid/yaspgrid.hh>
+
+typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+
+typedef testing::Types<Dune::Stuff::Functions::Spe10::Model2<DuneYaspGrid3dEntityType, double, 3, double, 3, 3>>
+ Spe10Model2FunctionYaspGridEntityTypes;
+
+TEST_STRUCT_GENERATOR(Spe10Model2FunctionTest, YaspGridEntity)
+TYPED_TEST_CASE(Spe10Model2FunctionTestYaspGridEntityTest, Spe10Model2FunctionYaspGridEntityTypes);
+TYPED_TEST(Spe10Model2FunctionTestYaspGridEntityTest, provides_required_methods)
+{
+ this->check();
+}
+
+//# if HAVE_ALUGRID
+//# include <dune/stuff/common/disable_warnings.hh>
+//# include <dune/grid/alugrid.hh>
+//# include <dune/stuff/common/reenable_warnings.hh>
+
+// typedef Dune::ALUGrid< 2, 2, Dune::simplex, Dune::nonconforming >::Codim< 0 >::Entity DuneAluSimplexGrid2dEntityType;
+// typedef Dune::ALUGrid< 3, 3, Dune::simplex, Dune::nonconforming >::Codim< 0 >::Entity DuneAluSimplexGrid3dEntityType;
+// typedef Dune::ALUGrid< 3, 3, Dune::cube, Dune::nonconforming >::Codim< 0 >::Entity DuneAluCubeGrid3dEntityType;
+
+// typedef testing::Types< Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 1,
+// 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 1, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 1, 3
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 2, 1
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 2, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 2, 3
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 3, 1
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 3, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid2dEntityType, double, 2, double, 3, 3
+// >
+
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 1, 1
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 1, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 1, 3
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 2, 1
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 2, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 2, 3
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 3, 1
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 3, 2
+// >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluSimplexGrid3dEntityType, double, 3, double, 3, 3
+// >
+
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 1, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 1, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 1, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 2, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 2, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 2, 3 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 3, 1 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 3, 2 >
+// , Dune::Stuff::Functions::Spe10::Model2< DuneAluCubeGrid3dEntityType, double, 3, double, 3, 3 >
+// > Spe10Model2FunctionAluGridEntityTypes;
+
+// TEST_STRUCT_GENERATOR(Spe10Model2FunctionTest, AluGridEntity)
+// TYPED_TEST_CASE(Spe10Model2FunctionTestAluGridEntityTest, Spe10Model2FunctionAluGridEntityTypes);
+// TYPED_TEST(Spe10Model2FunctionTestAluGridEntityTest, provides_required_methods) {
+// this->check();
+//}
+
+//# endif // HAVE_ALUGRID
+#endif // HAVE_DUNE_GRID