[function.spe10] added, not perfect yet

dune/stuff/function/spe10.hh

+#ifndef DUNE_STUFF_FUNCTION_SPE10_HH
+#define DUNE_STUFF_FUNCTION_SPE10_HH
+
+#include <sstream>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/stuff/common/parameter/tree.hh>
+#include <dune/stuff/common/color.hh>
+
+#include "interface.hh"
+
+
+namespace Dune {
+namespace Stuff {
+namespace Function {
+
+
+// forward, to allow for specialization
+template <class DomainFieldImp, int domainDim, class RangeFieldImp, int rangeDim>
+class Spe10Model1;
+
+
+template <class DomainFieldImp, class RangeFieldImp>
+class Spe10Model1<DomainFieldImp, 2, RangeFieldImp, 1> : public Interface<DomainFieldImp, 2, RangeFieldImp, 1>
+{
+public:
+  typedef Spe10Model1<DomainFieldImp, 2, RangeFieldImp, 1> ThisType;
+  typedef Interface<DomainFieldImp, 2, RangeFieldImp, 1> BaseType;
+
+  typedef typename BaseType::DomainFieldType DomainFieldType;
+  static const int dimDomain = BaseType::dimDomain;
+  typedef typename BaseType::DomainType DomainType;
+
+  typedef typename BaseType::RangeFieldType RangeFieldType;
+  static const int dimRange = BaseType::dimRange;
+  typedef typename BaseType::RangeType RangeType;
+
+  static std::string id()
+  {
+    return BaseType::id() + ".spe10.model1";
+  }
+
+private:
+  static const size_t numXelements = 100;
+  static const size_t numYelements = 1;
+  static const size_t numZelements = 20;
+
+public:
+  Spe10Model1(const std::string filename, const DomainType& _lowerLeft, const DomainType& _upperRight,
+              //              const std::vector< size_t >& _numElements,
+              const std::string _name = id(), const int _order = 0)
+    : lowerLeft_(_lowerLeft)
+    , upperRight_(_upperRight)
+    //    , numElements_(_numElements)
+    , name_(_name)
+    , order_(_order)
+  {
+    // sanity checks
+    std::stringstream msg;
+    size_t throw_up = 0;
+    for (int dd = 0; dd < dimDomain; ++dd) {
+      if (!(lowerLeft_[dd] < upperRight_[dd])) {
+        ++throw_up;
+        msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:") << " lowerLeft[" << dd << "] (" << lowerLeft_[dd]
+            << ") has to be smaller than upperRight[" << dd << "] (" << upperRight_[dd] << ")!";
+      }
+      //      if (!(numElements_[dd] > 0)) {
+      //        ++throw_up;
+      //        msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:")
+      //            << " numElements[" << dd << "] has to be positive (is " << numElements_[dd] << ")!";
+      //      }
+    } // for (int dd = 0; dd < dimDomain; ++dd)
+    //    if (!(numElements_[0] <= numXelements)) {
+    //      ++throw_up;
+    //      msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:")
+    //          << " numElements[0] has to be smaller than " << numXelements << " (is " << numElements_[0] << ")!";
+    //    }
+    //    if (!(numElements_[1] <= numZelements)) {
+    //      ++throw_up;
+    //      msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:")
+    //          << " numElements[1] has to be smaller than " << numZelements << " (is " << numElements_[1] << ")!";
+    //    }
+    // read all the data from the file
+    std::ifstream datafile(filename);
+    if (datafile.is_open()) {
+      static const size_t entriesPerDim = numXelements * numYelements * numZelements;
+      // create storage (there should be exactly 6000 values in the file)
+      data_           = new double[3 * entriesPerDim];
+      double tmpValue = 0.0;
+      size_t counter = 0;
+      while (datafile >> tmpValue && counter < 3 * entriesPerDim) {
+        data_[counter] = tmpValue;
+        ++counter;
+      }
+      datafile.close();
+      if (counter != 3 * entriesPerDim) {
+        ++throw_up;
+        msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:") << " wrong number of entries in '" << filename
+            << "' (are " << counter << ", should be " << 3 * entriesPerDim << ")!";
+      }
+    } else { // if (datafile)
+      ++throw_up;
+      msg << "\n" << Dune::Stuff::Common::colorStringRed("ERROR:") << " could not open '" << filename << "'!";
+    } // if (datafile)
+    // throw up, if needed
+    if (throw_up)
+      DUNE_THROW(Dune::RangeError, msg.str());
+    //    // keep the data we need
+    //    data_ = new double[numElements_[0] * numElements_[1]];
+
+
+    //    // clean up
+    //    delete data;
+  } // Spe10Model1
+
+  ~Spe10Model1()
+  {
+    delete data_;
+  }
+
+  static Dune::ParameterTree createSampleDescription(const std::string subName = "")
+  {
+    Dune::ParameterTree description;
+    description["filename"]   = "perm_case1.dat";
+    description["lowerLeft"]  = "[0.0; 0.0]";
+    description["upperRight"] = "[762.0; 15.24]";
+    //    description["numElements"] = "[100; 20]";
+    description["name"]  = id();
+    description["order"] = "0";
+    if (subName.empty())
+      return description;
+    else {
+      Dune::Stuff::Common::ExtendedParameterTree extendedDescription;
+      extendedDescription.add(description, subName);
+      return extendedDescription;
+    }
+  } // ... createSampleDescription(...)
+
+  static ThisType* createFromDescription(const Dune::Stuff::Common::ExtendedParameterTree description)
+  {
+    // get data
+    const std::string filenameIn                    = description.get<std::string>("filename");
+    const std::vector<DomainFieldType> lowerLeftIn  = description.getVector<DomainFieldType>("lowerLeft", dimDomain);
+    const std::vector<DomainFieldType> upperRightIn = description.getVector<DomainFieldType>("upperRight", dimDomain);
+    //    const std::vector< size_t >          numElements  = description.getVector< size_t >(         "numElements",
+    //                                                                                                 dimDomain);
+    const std::string nameIn = description.get<std::string>("name", id());
+    const int orderIn        = description.get<int>("order", 0);
+    // convert and leave the checks to the constructor
+    DomainType lowerLeft;
+    DomainType upperRight;
+    for (int dd = 0; dd < dimDomain; ++dd) {
+      lowerLeft[dd]  = lowerLeftIn[dd];
+      upperRight[dd] = upperRightIn[dd];
+    }
+    // create and return
+    return new ThisType(filenameIn, lowerLeft, upperRight, /*numElements,*/ nameIn, orderIn);
+  } // static ThisType createFromParamTree(const Dune::ParameterTree paramTree)
+
+  const DomainType& lowerLeft() const
+  {
+    return lowerLeft_;
+  }
+
+  const DomainType& upperRight() const
+  {
+    return upperRight_;
+  }
+
+  //  const std::vector< size_t >& numElements() const
+  //  {
+  //    return numElements_;
+  //  }
+
+  virtual std::string name() const
+  {
+    return name_;
+  }
+
+  virtual int order() const
+  {
+    return order_;
+  }
+
+  virtual void evaluate(const DomainType& x, RangeType& ret) const
+  {
+    // decide on the interval x belongs to
+    Dune::FieldVector<size_t, dimDomain> interval;
+    interval[0] = std::floor(numXelements * ((x[0] - lowerLeft_[0]) / (upperRight_[0] - lowerLeft_[0])));
+    interval[1] = std::floor(numZelements * ((x[1] - lowerLeft_[1]) / (upperRight_[1] - lowerLeft_[1])));
+    //    for (size_t dd = 0; dd < dimDomain; ++dd) {
+    //      interval[dd] = std::floor(numElements_[dd]*(           (x[dd] - lowerLeft_[dd])
+    //                                                  /(upperRight_[dd] - lowerLeft_[dd])));
+    //    }
+    //    if (interval[0] >= numElements_[0] || interval[1] >= numElements_[1]) {
+    if (interval[0] >= numXelements || interval[1] >= numZelements) {
+      ret[0] = 0;
+    } else {
+      const size_t index = interval[0] + numXelements * 0 + numXelements * numYelements * interval[1];
+      ret[0]             = data_[index];
+    }
+  } // virtual void evaluate(const DomainType& x, RangeType& ret) const
+
+private:
+  const std::string filename_;
+  const DomainType lowerLeft_;
+  const DomainType upperRight_;
+  //  const std::vector< size_t > numElements_;
+  const std::string name_;
+  const int order_;
+  double* data_;
+}; // class Spe10Model1< DomainFieldImp, 2, RangeFieldImp, 1 >
+
+
+} // namespace Function
+} // namespace Stuff
+} // namespace Dune
+
+#endif // DUNE_STUFF_FUNCTION_SPE10_HH