[!3] Resolve "Include more options to initializate stetes"

Merge branch '2-include-more-options-to-initializate-stetes' into 'master'

ref:copasi/dune-copasi

### What does this MR do?

Implements TIFF input data for the solver.

The section model.data is read if defined. It should contain a key assigned to
a path of a valid gayscale tiff of 16 bits. e.g.

    [model.data]
    tiff_data = data/tiff/flower-minisblack-16.tif

The key can be latter used in the initial condition of a compartment as a
function of the **grid** \$x\$ and \$y\$ coordinates. If arguments are out of
bound of the tiff image the function will return 0.

    [model.<compartment>.initial]
    u_0 = tiff_data(x,y)

The arguments of the function may be anything defined on muparser. This is
useful to shift and scale the data. Additionaly, data is interpreted between 0
and 1 depending on the photometric information (*zero_is_white* or
*zero_is_black*). In case data has another scale, its result may be scaled as
any other function in muparser. e.g.

    [model.<compartment>.initial]
    u_0 = 1.3+5.1*tiff_data(2*x, y+0.25)

### Is there something that needs to be double checked?

<!-- Is there something a reviewer should look out for _especially_? -->

No

### Can this MR be accepted?

-   [x] Implemented \...
    -   [x] Generic tiff reader
    -   [x] muParser function for tiff data
-   [x] Added/Updated tests:
    -   [x] Unit test for tiff reader
    -   [x] System test for initial condition
-   [ ] Pipelines passing <!-- please check for new warnings -->
     <!-- change all occurences of <branch> for your branch name -->

    -   [x] [![Build Status]]
    -   [x] [![Build Status][1]]
    -   [ ] [![Build status][2]]

-   [x] Delete branch option set <!-- unless there's a good reason -->

### Related issues

Closes [#2]

<!-- For automatic closing, do not forget the commas between issue numbers-->

<!--
PLEASE READ THIS!

A Merge Request should be associated to a certain task or issue.
Its changes are supposed to be merged into the master branch.

Briefly explain __how__ you achieved the proposal of the task.

IMPORTANT: Make sure to set the merge request WIP if you are not finished yet.
-->

See merge request [!3]

  [Build Status]: https://gitlab.dune-project.org/copasi/dune-copasi/badges/2-include-more-options-to-initializate-stetes/pipeline.svg
  [![Build Status]]: https://gitlab.dune-project.org/copasi/dune-copasi/pipelines
  [1]: https://travis-ci.org/SoilRos/dune-copasi.svg?branch=2-include-more-options-to-initializate-stetes
  [![Build Status][1]]: https://travis-ci.org/SoilRos/dune-copasi
  [2]: https://ci.appveyor.com/api/projects/status/6605joy2w17qvca8/branch/2-include-more-options-to-initializate-stetes?svg=true
  [![Build status][2]]: https://ci.appveyor.com/project/SoilRos/dune-copasi/history
  [#2]: gitlab.dune-project.org/NoneNone/issues/2
  [!3]: gitlab.dune-project.org/copasi/dune-copasi/merge_requests/3


Closes #2

.ci/README.md

-[![Build Status](https://gitlab.dune-project.org/santiago.ospina/dune-copasi/badges/master/pipeline.svg)](https://gitlab.dune-project.org/santiago.ospina/dune-copasi/pipelines)
+[![Build Status](https://gitlab.dune-project.org/copasi/dune-copasi/badges/master/pipeline.svg)](https://gitlab.dune-project.org/copasi/dune-copasi/pipelines)
 [![Build Status](https://travis-ci.org/SoilRos/dune-copasi.svg?branch=master)](https://travis-ci.org/SoilRos/dune-copasi)
 [![Build status](https://ci.appveyor.com/api/projects/status/6605joy2w17qvca8/branch/master?svg=true)](https://ci.appveyor.com/project/SoilRos/dune-copasi/branch/master)
 

.ci/setup.sh

@@ -82,7 +82,7 @@ cd dune-logging
 git apply ../dune-copasi/.ci/dune-logging.patch
 cd ../
 
-ls 
+ls
 
 # python virtual environment does not work in windows yet
 if [[ ! $MSYSTEM ]]; then

.travis.yml

@@ -32,7 +32,7 @@ before_install:
   - gcc --version
   - g++ --version
   - cd ..
-  # download muparser as static library
+  # download muparser, gmp and libtiff as static libraries
   - wget "https://github.com/lkeegan/libsbml-static/releases/latest/download/libsbml-static-$TRAVIS_OS_NAME.tgz"
   - tar xzvf libsbml-static-$TRAVIS_OS_NAME.tgz
 install:
@@ -42,17 +42,19 @@ install:
   - echo 'CMAKE_FLAGS+=" -DGMPXX_INCLUDE_DIR:PATH='"$PWD"'/gmp/include -DGMPXX_LIB:FILEPATH='"$PWD"'/gmp/lib/libgmpxx.a -DGMP_LIB:FILEPATH='"$PWD"'/gmp/lib/libgmp.a"' >> dune.opts
   - echo 'CMAKE_FLAGS+=" -DCMAKE_DISABLE_FIND_PACKAGE_QuadMath=TRUE -DBUILD_TESTING=OFF"' >> dune.opts
   - echo 'CMAKE_FLAGS+=" -DDUNE_USE_ONLY_STATIC_LIBS=ON -DCMAKE_BUILD_TYPE=Release"' >> dune.opts
-  - echo 'CMAKE_FLAGS+=" -DCMAKE_PREFIX_PATH='"$PWD"'/muparser -DF77=true "' >> dune.opts
-  - echo 'CMAKE_FLAGS+=" -DCMAKE_PREFIX_PATH='"$PWD"'/muparser -DF77=true "' >> dune.opts
+  - echo 'CMAKE_FLAGS+=" -DCMAKE_PREFIX_PATH='"$PWD"' -DF77=true "' >> dune.opts
+  - echo 'CMAKE_FLAGS+=" -Dmuparser_ROOT='"$PWD"'/muparser "' >> dune.opts
+  - echo 'CMAKE_FLAGS+=" -DTIFF_ROOT='"$PWD"'/libtiff "' >> dune.opts
   - if [ "$TRAVIS_OS_NAME" = "linux" ]; then echo 'CMAKE_FLAGS+="-DPYTHON_EXECUTABLE:FILEPATH=/usr/bin/python3"' >> dune.opts; fi
+  - echo 'MAKE_FLAGS="-j2 VERBOSE=1"' >> dune.opts
   - export DUNE_OPTIONS_FILE="dune.opts"
   - export DUNECONTROL=./dune-common/bin/dunecontrol
   - echo "${DUNE_OPTIONS_FILE}"
   # setup dependencies
   - bash dune-copasi/.ci/setup.sh
-script: 
+script:
   - bash dune-copasi/.ci/build.sh
-after_script: 
+after_script:
   - bash dune-copasi/.ci/unit_tests.sh
   - bash dune-copasi/.ci/system_tests.sh
   # check dependencies of executable

CMakeLists.txt

@@ -9,6 +9,8 @@ if(NOT (dune-common_DIR OR dune-common_ROOT OR
       ${PROJECT_BINARY_DIR})
 endif()
 
+cmake_policy(SET CMP0074 NEW)
+
 #find dune-common and set the module path
 find_package(dune-common REQUIRED)
 list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"

README.md

-[![Build Status](https://gitlab.dune-project.org/santiago.ospina/dune-copasi/badges/master/pipeline.svg)](https://gitlab.dune-project.org/santiago.ospina/dune-copasi/pipelines)
+[![Build Status](https://gitlab.dune-project.org/copasi/dune-copasi/badges/master/pipeline.svg)](https://gitlab.dune-project.org/copasi/dune-copasi/pipelines)
 [![Build Status](https://travis-ci.org/SoilRos/dune-copasi.svg?branch=master)](https://travis-ci.org/SoilRos/dune-copasi)
 [![Build status](https://ci.appveyor.com/api/projects/status/6605joy2w17qvca8/branch/master?svg=true)](https://ci.appveyor.com/project/SoilRos/dune-copasi/history)
 

appveyor.yml

@@ -6,15 +6,24 @@ clone_folder: C:\msys64\home\appveyor\dune-copasi
 
 install:
   - set PATH=C:\msys64\usr\bin;%PATH%
-  # download muparser & gmp as static libraries
+  # download muparser, gmp and libtiff as static libraries
   - appveyor DownloadFile "https://github.com/lkeegan/libsbml-static/releases/latest/download/libsbml-static-windows.zip"
   - 7z x libsbml-static-windows.zip
   - mv muparser C:\msys64\home\appveyor\.
   - mv gmp C:\msys64\home\appveyor\.
+  - mv libtiff C:\msys64\home\appveyor\.
 
 before_build:
   - set DEPS_DIR=/home/appveyor
-  - echo CMAKE_FLAGS="-G 'Unix Makefiles' -DCMAKE_VERBOSE_MAKEFILE=1 -DDUNE_PYTHON_VIRTUALENV_SETUP=0 -DDUNE_PYTHON_ALLOW_GET_PIP=0 -DDUNE_PYTHON_VIRTUALENV_PATH=%DEPS_DIR%/dune-python-venv  -DGMPXX_INCLUDE_DIR:PATH=%DEPS_DIR%/gmp/include -DGMPXX_LIB:FILEPATH=%DEPS_DIR%/gmp/lib/libgmpxx.a -DGMP_LIB:FILEPATH=%DEPS_DIR%/gmp/lib/libgmp.a -DCMAKE_DISABLE_FIND_PACKAGE_QuadMath=TRUE -DBUILD_TESTING=OFF -DDUNE_USE_ONLY_STATIC_LIBS=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH=%DEPS_DIR%/muparser -DF77=true -DCMAKE_CXX_FLAGS='-Wa,-mbig-obj -static -static-libgcc -static-libstdc++' " > C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS="-G 'Unix Makefiles' -DCMAKE_BUILD_TYPE=Release  -DCMAKE_VERBOSE_MAKEFILE=1 " > C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DDUNE_PYTHON_VIRTUALENV_SETUP=0 -DDUNE_PYTHON_ALLOW_GET_PIP=0" >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DDUNE_PYTHON_VIRTUALENV_PATH=%DEPS_DIR%/dune-python-venv" >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DGMPXX_INCLUDE_DIR:PATH=%DEPS_DIR%/gmp/include -DGMPXX_LIB:FILEPATH=%DEPS_DIR%/gmp/lib/libgmpxx.a" >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DGMP_LIB:FILEPATH=%DEPS_DIR%/gmp/lib/libgmp.a -DCMAKE_DISABLE_FIND_PACKAGE_QuadMath=TRUE" >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DBUILD_TESTING=OFF -DDUNE_USE_ONLY_STATIC_LIBS=ON" >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DCMAKE_PREFIX_PATH=%DEPS_DIR% -Dmuparser_ROOT=%DEPS_DIR%/muparser -DTIFF_ROOT=%DEPS_DIR%/libtiff " >> C:\msys64\home\appveyor\dune.txt
+  - echo CMAKE_FLAGS+=" -DF77=true -DCMAKE_CXX_FLAGS='-Wa,-mbig-obj -static -static-libgcc -static-libstdc++' " >> C:\msys64\home\appveyor\dune.txt
+  - echo MAKE_FLAGS="-j2 VERBOSE=1" >> C:\msys64\home\appveyor\dune.txt
   - type C:\msys64\home\appveyor\dune.txt
   - set DUNE_OPTIONS_FILE=dune.txt
   - set DUNECONTROL=./dune-common/bin/dunecontrol

cmake/modules/DuneCopasiMacros.cmake

 # File for module specific CMake tests.
 
-find_package (muparser REQUIRED)
\ No newline at end of file
+find_package (muparser REQUIRED)
+find_package (TIFF REQUIRED)

docker/dune-copasi.dockerfile

@@ -3,4 +3,4 @@ FROM ${BASE_IMAGE}
 ARG DUNECI_PARALLEL=2
 ARG BRANCH=master
 
-RUN duneci-install-module -b ${BRANCH} https://gitlab.dune-project.org/santiago.ospina/dune-copasi.git
\ No newline at end of file
+RUN duneci-install-module -b ${BRANCH} https://gitlab.dune-project.org/copasi/dune-copasi.git
\ No newline at end of file

dune/copasi/local_operator.hh

@@ -183,7 +183,7 @@ public:
           grid_view, d_eq);
       _reaction_gf[k] =
         std::make_shared<ExpressionToGridFunctionAdapter<GridView, RF>>(
-          grid_view, r_eq, reaction_keys);
+          grid_view, r_eq, true, reaction_keys);
 
       for (std::size_t l = 0; l < _lfs_components.size(); l++) {
         const auto j = _lfs_components.size() * k + l;
@@ -192,7 +192,7 @@ public:
 
         _jacobian_gf[j] =
           std::make_shared<ExpressionToGridFunctionAdapter<GridView, RF>>(
-            grid_view, j_eq, reaction_keys);
+            grid_view, j_eq, true, reaction_keys);
         if (k == l) {
           _component_pattern.insert(std::make_pair(k, l));
           continue;
@@ -238,6 +238,20 @@ public:
               pattern.addLink(lfsv.child(i), k, lfsu.child(j), l);
   }
 
+  void setTime(double t)
+  {
+    Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>::setTime(t);
+    for (auto& gf : _jacobian_gf)
+      if (gf)
+        gf->set_time(t);
+    for (auto& gf : _reaction_gf)
+      if (gf)
+        gf->set_time(t);
+    for (auto& gf : _diffusion_gf)
+      if (gf)
+        gf->set_time(t);
+  }
+
   template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
   void jacobian_apply_volume(const EG& eg,
                              const LFSU& lfsu,

dune/copasi/local_operator_multidomain.hh

@@ -235,6 +235,13 @@ public:
     }
   }
 
+  void setTime(double t)
+  {
+    Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>::setTime(t);
+    for (auto lp : _local_operator)
+      lp->setTime(t);
+  }
+
   template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
   void jacobian_apply_volume(const EG& eg,
                              const LFSU& lfsu,

dune/copasi/model_multidomain_diffusion_reaction.cc

@@ -11,6 +11,7 @@
 #include <dune/pdelab/gridfunctionspace/gridfunctionadapter.hh>
 
 #include <dune/copasi/model_multidomain_diffusion_reaction.hh>
+#include <dune/copasi/muparser_data_handler.hh>
 
 namespace Dune::Copasi {
 
@@ -38,6 +39,12 @@ ModelMultiDomainDiffusionReaction<Traits>::
   std::vector<std::shared_ptr<CompartementGridFunction>> comp_functions(
     _domains);
 
+  std::vector<std::shared_ptr<TIFFGrayscale<std::size_t>>> data(0);
+  std::vector<std::string> data_id(0);
+  MuParserDataHandler<TIFFGrayscale<unsigned short>> _mu_data_handler;
+  if (_config.hasSub("data"))
+    _mu_data_handler.add_tiff_functions(_config.sub("data"));
+
   for (auto& [op, state] : _states) {
     _logger.trace("interpolation of operator {}"_fmt, op);
     for (std::size_t i = 0; i < _domains; ++i) {
@@ -61,7 +68,11 @@ ModelMultiDomainDiffusionReaction<Traits>::
 
         _logger.trace("creating grid function for variable: {}"_fmt, var);
         std::string eq = intial_config[var];
-        functions[count] = std::make_shared<GridFunction>(_grid_view, eq);
+        functions[count] =
+          std::make_shared<GridFunction>(_grid_view, eq, false);
+        _mu_data_handler.set_functions(functions[count]->parser());
+        functions[count]->compile_parser();
+        functions[count]->set_time(current_time());
         count++;
       }
       // the second is because of the "ghost" child for empty compartments
@@ -237,7 +248,11 @@ ModelMultiDomainDiffusionReaction<Traits>::
     auto& lop = _local_operators[op];
     auto& tlop = _temporal_local_operators[op];
 
-    MBE mbe((int)pow(3, dim));
+    std::size_t max_comps(0);
+    for (std::size_t i = 0; i < gfs->degree(); i++)
+      max_comps = std::max(max_comps,gfs->child(i).degree());
+    
+    MBE mbe((int)pow(3, dim)*max_comps);
 
     _logger.trace("create spatial grid operator {}"_fmt, op);
     _spatial_grid_operators[op] = std::make_shared<GOS>(
@@ -371,7 +386,10 @@ ModelMultiDomainDiffusionReaction<Traits>::step()
   double max_error = std::numeric_limits<double>::max();
   auto states_after = _states;
 
-  _logger.info("Time Step {:.2e} + {:.2e} -> {:.2e}"_fmt, current_time(), dt, current_time()+dt);
+  _logger.info("Time Step {:.2e} + {:.2e} -> {:.2e}"_fmt,
+               current_time(),
+               dt,
+               current_time() + dt);
   std::size_t op_iter = 0;
   do {
     const auto states_before = states_after;

dune/copasi/muparser_data_handler.hh

+#ifndef DUNE_MUPARSER_DATA_HANDLER_HH
+#define DUNE_MUPARSER_DATA_HANDLER_HH
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/parametertree.hh>
+
+#include <muParser.h>
+
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace Dune::Copasi {
+
+// This is a trick to make happy muparser about static functions
+
+template<class T>
+struct MuParserDataHandler
+{
+  MuParserDataHandler() {}
+
+  ~MuParserDataHandler()
+  {
+    MuParserDataHandler<T>::_names.clear();
+    MuParserDataHandler<T>::_functions.clear();
+  }
+
+  template<int i>
+  static double function_wrapper(double x, double y)
+  {
+    assert(MuParserDataHandler<T>::_functions[i]);
+    return (*MuParserDataHandler<T>::_functions[i])(x, y);
+  }
+
+  void add_tiff_functions(const ParameterTree& data_config)
+  {
+    const auto& keys = data_config.getValueKeys();
+    for (std::size_t i = 0; i < keys.size(); i++) {
+      _names.push_back(keys[i]);
+      std::string filename = data_config[keys[i]];
+      _functions.push_back(std::make_shared<T>(filename));
+    }
+  }
+
+  template<std::size_t max_data = 100>
+  void set_functions(mu::Parser& parser)
+  {
+    auto indices = Dune::range(std::integral_constant<std::size_t, max_data>{});
+    Dune::Hybrid::forEach(indices, [&](auto i) {
+      if (i < _functions.size())
+        parser.DefineFun(_names[i], function_wrapper<i>);
+    });
+  }
+
+  static std::vector<std::string> _names;
+  static std::vector<std::shared_ptr<T>> _functions;
+};
+
+template<class T>
+std::vector<std::string> MuParserDataHandler<T>::_names = {};
+
+template<class T>
+std::vector<std::shared_ptr<T>> MuParserDataHandler<T>::_functions = {};
+
+} // namespace Dune::Copasi
+
+#endif // DUNE_MUPARSER_DATA_HANDLER_HH
\ No newline at end of file

dune/copasi/pdelab_expression_adapter.hh

 #ifndef DUNE_COPASI_GRID_FUNCTION_EXPRESSION_ADAPTER_HH
 #define DUNE_COPASI_GRID_FUNCTION_EXPRESSION_ADAPTER_HH
 
+#include <dune/copasi/tiff_grayscale.hh>
+
 #include <dune/pdelab/common/function.hh>
 
 #include <dune/logging/logging.hh>
@@ -12,6 +14,7 @@
 
 #include <algorithm>
 #include <string>
+#include <type_traits>
 
 namespace Dune::Copasi {
 
@@ -25,12 +28,17 @@ public:
   using Traits = PDELab::GridFunctionTraits<GV, RF, 1, FieldVector<RF, 1>>;
 
   //! construct from grid view
+
   ExpressionToGridFunctionAdapter(const GV& grid_view,
                                   const std::string& equation,
+                                  bool do_compile_parser = true,
                                   std::vector<std::string> other_variables = {})
     : _logger(Logging::Logging::componentLogger({}, "model"))
     , _gv(grid_view)
+    , _time(0.)
     , _other_value(other_variables.size())
+    , _expr(equation)
+    , _compiled(false)
   {
 
     constexpr int dim = Traits::dimDomain;
@@ -41,25 +49,21 @@ public:
     _parser.DefineConst("pi", StandardMathematicalConstants<double>::pi());
     _parser.DefineConst("dim", dim);
 
+    _parser.DefineVar("t", &_time);
     _parser.DefineVar("x", &_pos_global[0]);
     _parser.DefineVar("y", &_pos_global[1]);
+
     if constexpr (dim == 3)
       _parser.DefineVar("z", &_pos_global[2]);
 
-    // set up parser expression
-    try {
-      for (size_t i = 0; i < other_variables.size(); i++) {
-        _logger.trace("define extra variable: {}"_fmt, other_variables[i]);
-        _parser.DefineVar(other_variables[i], &_other_value[i]);
-      }
-      _logger.trace("compile expression: {}"_fmt, equation);
-      _parser.SetExpr(equation);
-      // try to evaluate once
-      _parser.Eval();
-    } catch (mu::Parser::exception_type& e) {
-      handle_parser_error(e);
+    for (size_t i = 0; i < other_variables.size(); i++) {
+      _logger.trace("define extra variable: {}"_fmt, other_variables[i]);
+      _parser.DefineVar(other_variables[i], &_other_value[i]);
     }
 
+    if (do_compile_parser)
+      compile_parser();
+
     _logger.debug("ExpressionToGridFunctionAdapter constructed"_fmt);
   }
 
@@ -75,6 +79,7 @@ public:
   template<class E, class D, class R>
   void evaluate(const E& e, const D& x, R& y) const
   {
+    assert(_compiled);
     // update position storage
     _pos_global = e.geometry().global(x);
 
@@ -92,6 +97,27 @@ public:
     _other_value = other_value;
   }
 
+  void compile_parser()
+  {
+    try {
+      _logger.trace("compile expression: {}"_fmt, _expr);
+      _parser.SetExpr(_expr);
+      // try to evaluate once
+      _parser.Eval();
+    } catch (mu::Parser::exception_type& e) {
+      handle_parser_error(e);
+    }
+    _compiled = true;
+  }
+
+  mu::Parser& parser()
+  {
+    assert(not _compiled);
+    return _parser;
+  }
+
+  void set_time(double t) { _time = t; }
+
 private:
   /// Output information on the parser error and throw DUNE exception
   /**
@@ -116,10 +142,15 @@ private:
   /// Cache for the evaluation position
   mutable typename Traits::DomainType _pos_global;
 
+  double _time;
+
   DynamicVector<RF> _other_value;
 
   /// The parser instance
   mu::Parser _parser;
+
+  std::string _expr;
+  bool _compiled;
 };
 
 } // namespace Dune::Copasi

dune/copasi/tiff_grayscale.hh

+#ifndef DUNE_COPASI_TIFF_GRAYSCALE_HH
+#define DUNE_COPASI_TIFF_GRAYSCALE_HH
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+
+#include <tiffio.h>
+
+#include <string>
+#include <queue>
+#include <memory>
+
+namespace Dune::Copasi {
+
+template<class T>
+class TIFFGrayscale
+{
+  static_assert(std::is_integral_v<T>, "T must be an integral type");
+  static_assert(std::is_unsigned_v<T>, "T must be an unsigned type");
+
+  struct TIFFGrayscaleRow
+  {
+    TIFFGrayscaleRow(TIFF* const tiff_file, const T& row, const short& col_size, const bool& zero)
+      : _row(row)
+      , _col_size(col_size)
+      , _zero(zero)
+    {
+      T* raw_buffer = (T*)_TIFFmalloc(TIFFScanlineSize(tiff_file));
+      auto deleter = [](auto& ptr){_TIFFfree(ptr);};
+      _tiff_buffer = std::shared_ptr<T>(raw_buffer, deleter);
+      TIFFReadScanline(tiff_file, _tiff_buffer.get(), _row);
+    }
+
+    double operator[](const T& col) const
+    {
+      assert((short)col < _col_size);
+      const T max = std::numeric_limits<T>::max();
+      T val = *(_tiff_buffer.get() + col);
+      val = _zero ? val : max-val; 
+      return (double)val/max;
+    }
+
+    std::size_t size() const { return static_cast<std::size_t>(_col_size); }
+    std::size_t row() const { return static_cast<std::size_t>(_row); }
+
+  private:
+    std::shared_ptr<T> _tiff_buffer;
+    const T _row;
+    const short _col_size;
+    const bool _zero;
+  };
+
+public:
+  TIFFGrayscale(const std::string& filename)
+    : _tiff_file(TIFFOpen(filename.c_str(), "r"))
+  {
+    if (not _tiff_file)
+      DUNE_THROW(IOError, "Error opening TIFF file '" << filename << "'.");
+
+    short photometric;
+    TIFFGetField(_tiff_file, TIFFTAG_PHOTOMETRIC, &photometric);
+    if ((photometric != PHOTOMETRIC_MINISWHITE) and
+        (photometric != PHOTOMETRIC_MINISBLACK))
+      DUNE_THROW(IOError,
+                 "TIFF file '" << filename << "' must be in grayscale.");
+    _zero = (bool)photometric;
+
+    short bits_per_sample;
+    TIFFGetField(_tiff_file, TIFFTAG_BITSPERSAMPLE, &bits_per_sample);
+
+    if (bits_per_sample != 8 * sizeof(T)) {
+      TIFFClose(_tiff_file);
+      DUNE_THROW(IOError,
+                 "TIFF file '" << filename
+                               << "' contains a non-readable grayscale field.");
+    }
+
+    TIFFGetField(_tiff_file, TIFFTAG_IMAGELENGTH, &_row_size);
+    TIFFGetField(_tiff_file, TIFFTAG_IMAGEWIDTH, &_col_size);
+    TIFFGetField(_tiff_file, TIFFTAG_XRESOLUTION, &_x_res);
+    TIFFGetField(_tiff_file, TIFFTAG_YRESOLUTION, &_y_res);
+    assert(FloatCmp::gt(_x_res, 0.f));
+    assert(FloatCmp::gt(_y_res, 0.f));
+    _x_off = _y_off = 0.;
+    TIFFGetField(_tiff_file, TIFFTAG_XPOSITION, &_x_off);
+    TIFFGetField(_tiff_file, TIFFTAG_YPOSITION, &_y_off);
+  }
+
+  ~TIFFGrayscale()
+  {
+    TIFFClose(_tiff_file);
+  }
+
+  // if rows are read concurrently, this object has to be copied
+  const TIFFGrayscaleRow& operator[](T row) const
+  {
+    assert((short)row < _row_size);
+    return cache(row);
+  }
+
+  template<class DF>
+  double operator()(const DF& x, const DF& y)
+  {
+    // here we assume TIFFTAG_RESOLUTIONUNIT has same units as the grid.
+    // since we never check grid units, we also do not check tiff units
+
+    // return 0 if not in the domain
+    if (FloatCmp::lt((float)x, _x_off) or FloatCmp::lt((float)y, _y_off)) 
+      return 0;
+
+    const T i = _x_res * (_x_off + x);
+    const T j = _row_size - _y_res * (_y_off + y) - 1;
+
+    // return 0 if not in the domain
+    if (i>=cols() or j>=rows()) 
+      return 0;
+    else
+      return (*this)[j][i];
+  }
+
+  std::size_t size() const { return cols(); }
+
+  std::size_t rows() const { return static_cast<std::size_t>(_row_size); }
+
+  std::size_t cols() const { return static_cast<std::size_t>(_col_size); }
+
+private:
+
+const TIFFGrayscaleRow& cache(T row) const
+{
+  auto it = _row_cache.rbegin();
+  while ((it != _row_cache.rend()) and (it->row() != row))
+    it++;
+  
+  if (it != _row_cache.rend())
+    return *it;
+  else
+    _row_cache.emplace_back(_tiff_file, row, _col_size, _zero);
+  
+  if (_row_cache.size() >= 8)
+    _row_cache.pop_front();
+
+  return *(_row_cache.rbegin());
+}
+
+private:
+  TIFF* _tiff_file;
+  mutable std::deque<TIFFGrayscaleRow> _row_cache;
+  short _row_size;
+  short _col_size;
+  float _x_res, _x_off;
+  float _y_res, _y_off;
+  bool _zero;
+};
+
+} // namespace Dune::Copasi
+
+#endif // DUNE_COPASI_TIFF_GRAYSCALE_HH
\ No newline at end of file

lib/CMakeLists.txt

 add_library(dune_copasi_lib STATIC ../dune/copasi/model_base.cc
                                    model_multidomain_diffusion_reaction.cc)
 
-target_link_libraries(dune_copasi_lib PUBLIC ${DUNE_LIBS} muparser::muparser)
+target_link_libraries(dune_copasi_lib PUBLIC ${DUNE_LIBS} TIFF::TIFF muparser::muparser)
 
 install(TARGETS dune_copasi_lib 
         ARCHIVE

src/CMakeLists.txt

@@ -2,4 +2,4 @@ add_executable(dune_copasi dune_copasi.cc)
 target_link_libraries(dune_copasi PRIVATE dune_copasi_lib)
 
 install(TARGETS dune_copasi RUNTIME
-        DESTINATION "${CMAKE_INSTALL_BINDIR}/")
\ No newline at end of file
+        DESTINATION "${CMAKE_INSTALL_BINDIR}/")

test/CMakeLists.txt

@@ -2,6 +2,7 @@ enable_testing()
 
 dune_symlink_to_source_files(FILES "grids" DESTINATION "." )
 dune_symlink_to_source_files(FILES "reference" DESTINATION "." )
+dune_symlink_to_source_files(FILES "data" DESTINATION "." )
 
 dune_add_test(NAME            test_dynamic_local_finite_element
               SOURCES         test_dynamic_local_finite_element.cc
@@ -19,6 +20,11 @@ dune_add_test(NAME            test_concepts_typetree
               SOURCES         test_concepts_typetree.cc
               LABELS          unit)
 
+dune_add_test(NAME            test_tiff_grayscale
+              SOURCES         test_tiff_grayscale.cc
+              LINK_LIBRARIES  TIFF::TIFF
+              LABELS          unit)
+
 if (dune-testtools_FOUND)
 
   add_system_test_per_target(
@@ -31,6 +37,11 @@ if (dune-testtools_FOUND)
     SCRIPT          dune_vtkcompare.py
     INIFILE         test_exp.mini)
 
+    add_system_test_per_target(
+      TARGET          dune_copasi
+      SCRIPT          dune_vtkcompare.py
+      INIFILE         test_tiff.mini)
+
   dune_add_system_test(
     SOURCE          test_jacobian.cc
     BASENAME        dune_copasi_test_jacobian