diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b673e1c0e796ae612baeb3832611d17fb3755f62..37a4f997b9446fc0a2f6f4263f23d429f49ef969 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -125,7 +125,6 @@ stages:
only:
- master
-
# debian gcc
setup:debian_gcc:
<<: *debian_gcc
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7113ee31e3ab8dcdd05f2a839341c9c524bf033a..4ecfbae737589684010b03ae64bfa9dcd470db18 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -26,6 +26,9 @@ Types of changes
-->
## [Unreleased]
+### Added
+- Code documentation
+- Brief installation instructions
## [0.1.0] - 2019-10-11
### Added
diff --git a/README.md b/README.md
index e76e7106d744b80aeb1dc8bc14ce77180a41da62..db6258af267b33c3cf0a3c03b87f0e519cc5be1d 100644
--- a/README.md
+++ b/README.md
@@ -2,88 +2,230 @@
[](https://travis-ci.org/SoilRos/dune-copasi)
[](https://ci.appveyor.com/project/SoilRos/dune-copasi/history)
-#### Dependencies
+# dune-copasi
+
+Solver for reaction-diffusion systems in multiple compartments
+
+ * Solve a reaction-diffusion system for each comartment
+ * Each compartment may have different system with different number of variables
+ * Neumann flux at the interface of compartments for variables with
+ the same name on the two compartments
+ * Easy to modify configuration file
+ * Initial conditions might be either TIFF file or a math expression.
+ * Solved using the finite element method
+ * Output in the VTK format
+ * Currently it only supports 2D simulations
+
+This project is made under the umbrella of the
+[*Distributed and Unified Numerics Environment* `DUNE`](https://www.dune-project.org/) and the
+[*Biochemical System Simulator* `COPASI`](http://copasi.org/).
+Altought the rationale of the design is always driven by biochemical process (e.g. cell biology),
+this software is not limited to this scope and can be used for other processes involving reaction-diffusion systems.
+
+## Graphical User Interface for SMBL files
+
+For those working in bio-informatics there exist a grafical user interface for SMBL files!
+The GUI is able to convert non-spatial SBML models of bio-chemical reactions into
+2d spatial models, and to simulate them with `dune-copasi`.
+
+https://github.com/lkeegan/spatial-model-editor
+
+## Installation
+
+This requires that you have installed the following packages before the actual installation of `dune-copasi`.
| Software | Version/Branch | Comments |
| ---------| -------------- | -------- |
-| muParser | - |
-| CMake | 3.10.2 |
+| [CMake](https://cmake.org/) | 3.1 |
| C++ compiler | [C++17](https://en.wikipedia.org/wiki/List_of_compilers#C++_compilers) |
-| [dune-common](https://gitlab.dune-project.org/santiago.ospina/dune-common) | support/dune-copasi
-| [dune-logging](https://gitlab.dune-project.org/staging/dune-logging) | master
-| [dune-geometry](https://gitlab.dune-project.org/core/dune-geometry) | master
-| [dune-grid](https://gitlab.dune-project.org/core/dune-grid) | master
-| [dune-uggrid](https://gitlab.dune-project.org/staging/dune-uggrid) | master
-| [dune-istl](https://gitlab.dune-project.org/core/dune-istl) | master
-| [dune-localfunctions](https://gitlab.dune-project.org/core/dune-localfunctions) | master
-| [dune-functions](https://gitlab.dune-project.org/staging/dune-functions) | master
-| [dune-typetree](https://gitlab.dune-project.org/santiago.ospina/dune-typetree) | support/dune-copasi
-| [dune-pdelab](https://gitlab.dune-project.org/santiago.ospina/dune-pdelab) | support/dune-copasi
-| [dune-multidomaingrid](https://gitlab.dune-project.org/santiago.ospina/dune-multidomaingrid) | support/dune-copasi
-
-<!--
-Preparing the Sources
-=========================
-
-Additional to the software mentioned in README you'll need the
-following programs installed on your system:
+| [libTIFF](http://www.libtiff.org/) | 3.6.1 |
+| [muParser](https://beltoforion.de/article.php?a=muparser) | 2.2.5 |
+| [dune-common](https://gitlab.dune-project.org/santiago.ospina/dune-common) | support/dune-copasi | https://gitlab.dune-project.org/santiago.ospina/dune-common
+| [dune-logging](https://gitlab.dune-project.org/staging/dune-logging) | master (recursive) | https://gitlab.dune-project.org/staging/dune-logging
+| [dune-geometry](https://gitlab.dune-project.org/core/dune-geometry) | master | https://gitlab.dune-project.org/core/dune-geometry
+| [dune-grid](https://gitlab.dune-project.org/core/dune-grid) | master | https://gitlab.dune-project.org/core/dune-grid
+| [dune-uggrid](https://gitlab.dune-project.org/staging/dune-uggrid) | master | https://gitlab.dune-project.org/staging/dune-uggrid
+| [dune-istl](https://gitlab.dune-project.org/core/dune-istl) | master | https://gitlab.dune-project.org/core/dune-istl
+| [dune-localfunctions](https://gitlab.dune-project.org/core/dune-localfunctions) | master | https://gitlab.dune-project.org/core/dune-localfunctions
+| [dune-functions](https://gitlab.dune-project.org/staging/dune-functions) | master | https://gitlab.dune-project.org/staging/dune-functions
+| [dune-typetree](https://gitlab.dune-project.org/santiago.ospina/dune-typetree) | support/dune-copasi | https://gitlab.dune-project.org/santiago.ospina/dune-typetree
+| [dune-pdelab](https://gitlab.dune-project.org/santiago.ospina/dune-pdelab) | support/dune-copasi | https://gitlab.dune-project.org/santiago.ospina/dune-pdelab
+| [dune-multidomaingrid](https://gitlab.dune-project.org/santiago.ospina/dune-multidomaingrid) | support/dune-copasi | https://gitlab.dune-project.org/santiago.ospina/dune-multidomaingrid
+
+The first four can be obtained by your prefered package manager in unix-like operating systems. e.g.
-```
- cmake >= 2.8.12
+```bash
+# if you are in a debian/ubuntu OS
+apt update
+apt install cmake gcc g++ libtiff-dev libmuparser-dev git
+
+# if you are in a macOS
+xcode-select --install # Apple Command Line Tools
+brew update
+brew install cmake gcc libtiff muparser git
```
-Getting started
----------------
+Now, the dune modules (including `dune-copasi`) can be all checkout in a same folder and be installed in one go.
-If these preliminaries are met, you should run
-
-```
- dunecontrol all
+```bash
+# prepare a folder to download and build dune modules
+mkdir ~/dune-modules && cd ~/dune-modules
+
+# fetch dependencies & dune-copasi in ~/dune-modules folder
+git clone -b support/dune-copasi https://gitlab.dune-project.org/santiago.ospina/dune-common
+git clone -b master --recursive https://gitlab.dune-project.org/staging/dune-logging
+git clone -b master https://gitlab.dune-project.org/core/dune-geometry
+git clone -b master https://gitlab.dune-project.org/core/dune-grid
+git clone -b master https://gitlab.dune-project.org/staging/dune-uggrid
+git clone -b master https://gitlab.dune-project.org/core/dune-istl
+git clone -b master https://gitlab.dune-project.org/core/dune-localfunctions
+git clone -b master https://gitlab.dune-project.org/staging/dune-functions
+git clone -b support/dune-copasi https://gitlab.dune-project.org/santiago.ospina/dune-typetree
+git clone -b support/dune-copasi https://gitlab.dune-project.org/santiago.ospina/dune-pdelab
+git clone -b support/dune-copasi https://gitlab.dune-project.org/santiago.ospina/dune-multidomaingrid
+git clone -b master https://gitlab.dune-project.org/copasi/dune-copasi
+
+# configure and build dune modules
+./dune-common/bin/dunecontrol make all
+
+# install dune-copasi (this operation may requiere sudo)
+./dune-common/bin/dunecontrol --only=dune-copasi bexec make install
+
+# if you do not want to install dune-copasi system-wide, you can set
+# the CMAKE_INSTALL_PREFIX to a non restricted folder
+# see https://cmake.org/cmake/help/latest/variable/CMAKE_INSTALL_PREFIX.html
+
+# remove source and build files
+cd .. && rm -r ~/dune-modules
```
-which will find all installed dune modules as well as all dune modules
-(not installed) which sources reside in a subdirectory of the current
-directory. Note that if dune is not installed properly you will either
-have to add the directory where the `dunecontrol` script resides (probably
-`./dune-common/bin`) to your path or specify the relative path of the script.
+For further info on dune module installation process, please check out
+the [dune-project web page](https://www.dune-project.org/doc/installation/)
-Most probably you'll have to provide additional information to `dunecontrol`
-(e. g. compilers, configure options) and/or make options.
+## Usage
-The most convenient way is to use options files in this case. The files
-define four variables:
+If you installed `dune-copasi` system-wide, you should be able to call the program
+`dune_copasi` from your command line accompained with a configuration file.
-```
-CMAKE_FLAGS flags passed to cmake (during configure)
+```bash
+dune_copasi config.ini
```
-An example options file might look like this:
+### Configuration File
-```bash
-#use this options to configure and make if no other options are given
-CMAKE_FLAGS=" \
--DCMAKE_CXX_COMPILER=g++-5 \
--DCMAKE_CXX_FLAGS='-Wall -pedantic' \
--DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
-```
+The configuration file follows [INI file format](https://en.wikipedia.org/wiki/INI_file).
+It should contain at least two sections: `grid`, and `model`, whereas a third section
+`logging` is optional.
-If you save this information into example.opts you can pass the opts file to
-dunecontrol via the `--opts option`, e. g.
+#### Grid
-```bash
- dunecontrol --opts=example.opts all
+The grid section is fairly simple as it only contains the path to a [gmsh file](http://gmsh.info/)
+and its initial refinement level:
+
+```ini
+[grid]
+file = my_gmsh_file.msh
+initial_level = 1
```
-More info
----------
+#### Model
-See
+The model section starts with the definitions of the simulation time interval
+and the polynomal order of the finite element method (currently only supports `1` and `2`).
-```bash
- dunecontrol --help
+```ini
+[model]
+begin_time = 0.
+end_time = 10
+time_step = 0.1
+order = 1
+```
+
+The following is the definition of the compartments of the model.
+Each compartment corresponds to a *physical group* in the gmsh jargon.
+Although the gmsh format allows you to name such physical groups,
+we still need to assign them to a `dune-copasi` compartmet and for that
+we use the *physical group* index. Notice that uses 0-based indices while
+`gmsh` uses 1-based indices. In other words,
+`<gmsh_physical_group> = <dune_copasi_compartment> - 1`.
+Let's say for example that there is two *physical groups* in our gmsh file
+and we are going to use them as `nucleous` and `cytoplasm` compartments:
+
+```ini
+[model.compartments]
+ # nucleous corresponds to the physical group 1 in the gmsh file
+nucleous = 0
+ # cytoplasm corresponds to the physical group 2 in the gmsh file
+cytoplasm = 0
```
-for further options.
+Now, each of these compartments will define its own initial conditions,
+its diffusion-reaction system, and its vtk writter. For that, you have to expand
+the `model` section with the defined compartments, e.g. `model.nucleous` or `model.cytoplasm`.
+The subsection `initial`, `reaction`, `diffusion` and `operator` define the system variables
+and its properties. You can put as many variables as desired as long as they are the same
+in this three subsections. Notice that each variable defines a new diffusion-reaction partial
+differential equation associated with it.
+
+ * The `initial` subsection allow the initialization of each of the variables.
+ * The `diffusion` subsection defines the diffusion coefficient math
+ expression associated with each variable. It may only depend on the grid coordinates `x` and `y`.
+ * The `reaction` subsection defines the reaction part of each equation in the PDE.
+ Since this is the souce of non-linearities, it allows to be dependent on other defined variables
+ within the compartment. This section has to include yet another subsection called `jacobian`.
+ * The `reaction.jacobian` subsection must describe the
+ [jacobian matrix](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant)
+ of the `reaction` part. It must follow the syntax of `d<var_i>_d<var_j>`, which
+ reads as the *partial derivate of `<var_i>` with respect to `<var_j>`*.
+ * The `operator` subsection is an experimental feature and we recommend to set all variables to the
+ same index, e.g. 0.
+ * Finally, the subsection `writer` will define the file name for the vtk output.
+
+For example, the following `mode.nucleous` section defines a [Gray-Scott
+model with `F=0.0420` and `k=0.0610`](http://mrob.com/pub/comp/xmorphia/F420/F420-k610.html):
+
+```ini
+[model.nucleous.initial]
+u = 0.5
+v = (x>0) && (x<0.5) && (y>0.) && (y<0.5) ? 1 : 0
+
+[model.nucleous.diffusion]
+u = 2e-5
+v = 2e-5/2
+
+[model.nucleous.reaction]
+u = -u*v^2+0.0420*(1-u)
+v = u*v^2-(0.0420+0.0610)*v
+
+[model.nucleous.reaction.jacobian]
+du_du = -v^2-0.0420
+du_dv = -2*u*v
+dv_du = v^2
+dv_dv = 2*u*v-(0.0420+0.0610)
+
+[model.nucleous.operator]
+u = 0
+v = 0
+
+[model.nucleous.writer]
+file_name = nucleous_output
+```
+The `model.cytoplasm` would have to be defined in similar way. An important aspect when working
+with different compartments is the interface fluxes. In `dune-copasi` thex fluxes are set
+automatically to [dirichlet-dirichlet](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition)
+boundary conditions iff the variable is shared between the two intersecting domains. Further
+improvements will come for interface treatment.
+
+#### Logging
+
+The logging settings are directly forwarded to the `dune-logging` module. Please check
+its doxygen documentation for detailed information. A simple configuration is the following:
+```ini
+[logging]
+# possible levels: off, critical, error, waring, notice, info, debug, trace, all
+default.level = trace
-The full build system is described in the `dune-common/doc/buildsystem` (Git version) or under `share/doc/dune-common/buildsystem` if you installed DUNE! -->
\ No newline at end of file
+[logging.sinks.stdout]
+pattern = [{reldays:0>2}-{reltime:8%T}] [{backend}] {msg}
+```
\ No newline at end of file
diff --git a/dune/copasi/CMakeLists.txt b/dune/copasi/CMakeLists.txt
index 37f27f3e2103ab36e8b3c7856f637fb82e6728b2..3dbb4df0941a249320ed12876df55db539ac4611 100644
--- a/dune/copasi/CMakeLists.txt
+++ b/dune/copasi/CMakeLists.txt
@@ -1,25 +1,5 @@
+add_subdirectory(common)
add_subdirectory(concepts)
-
-install(FILES coefficient_mapper.hh
- dynamic_local_basis.hh
- dynamic_local_coefficients.hh
- dynamic_local_finite_element.hh
- dynamic_local_interpolation.hh
- enum.hh
- gmsh_reader.hh
- grid_function_writer.hh
- local_operator_multidomain.hh
- local_operator.hh
- model_base.hh
- model_diffusion_reaction.cc
- model_diffusion_reaction.hh
- model_multidomain_diffusion_reaction.cc
- model_multidomain_diffusion_reaction.hh
- model_state.hh
- multidomain_entity_transformation.hh
- multidomain_local_finite_element_map.hh
- muparser_data_handler.hh
- pdelab_callable_adapter.hh
- pdelab_expression_adapter.hh
- tiff_grayscale.hh
- DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/dune/copasi")
+add_subdirectory(finite_element)
+add_subdirectory(grid)
+add_subdirectory(model)
diff --git a/dune/copasi/common/CMakeLists.txt b/dune/copasi/common/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0f8d8e5f28ada1bfa23fea434caa7f46eb9b26ef
--- /dev/null
+++ b/dune/copasi/common/CMakeLists.txt
@@ -0,0 +1,6 @@
+install(FILES coefficient_mapper.hh
+ enum.hh
+ muparser_data_handler.hh
+ pdelab_expression_adapter.hh
+ tiff_grayscale.hh
+ DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/dune/copasi/common")
\ No newline at end of file
diff --git a/dune/copasi/coefficient_mapper.hh b/dune/copasi/common/coefficient_mapper.hh
similarity index 52%
rename from dune/copasi/coefficient_mapper.hh
rename to dune/copasi/common/coefficient_mapper.hh
index 260dd832d3156b03c9eca33cc5f8d85f98cfffa6..4f95bddf78080a33f8773b4e6a298e03d0997c75 100644
--- a/dune/copasi/coefficient_mapper.hh
+++ b/dune/copasi/common/coefficient_mapper.hh
@@ -14,12 +14,38 @@
namespace Dune::Copasi {
+/**
+ * @brief Default coefficinet mapper
+ * @details Does nothing but evaluate the second an the third arguments of
+ * the parenthesis call into the first one. Useful to keep code
+ * compatible with other coefficient mappers
+ */
struct DefaultCoefficientMapper
{
+ /**
+ * @brief Constructor. Arguments in this class do nothimg.
+ *
+ * @param[in] operator_config The operator configuration file
+ * @param[in] this_op The operator for local coefficients
+ */
DefaultCoefficientMapper(const ParameterTree& operator_config,
std::size_t this_op)
{}
+ /**
+ * @brief Evaluates the coefficents for a give component and a degree of
+ * freedom. Either from a local vector x_view_local, or by other
+ * stored coefficients vector.
+ *
+ * @param[in] x_view_local The view on the local coefficients vector
+ * @param[in] comp The component
+ * @param[in] dof The degree of freedom
+ *
+ * @tparam XViewLocal The type fo the local coefficient view
+ *
+ * @return The coefficient value for the requested component and the
+ * degree of freedom
+ */
template<class XViewLocal>
auto operator()(const XViewLocal& x_view_local,
const std::size_t& comp,
@@ -28,17 +54,41 @@ struct DefaultCoefficientMapper
return x_view_local(comp, dof);
}
+ /**
+ * @brief Bind internal coefficients to an entity
+ *
+ * @param[in] <unnamed> The entity
+ *
+ * @tparam E The entity
+ */
template<class E>
void bind(const E&)
{}
+ /**
+ * @brief Unbind this coefficients to an entity
+ */
void unbind() {}
+ /**
+ * @brief Updates the internal coefficients with a exteral source.
+ *
+ * @param[in] <unnamed> Map to the model states
+ *
+ * @tparam T Map to the model states
+ */
template<class T>
void update(const T&)
{}
};
+/**
+ * @brief Model coefficinet mapper base
+ * @details Keeps a reference to external coefficients so that they are
+ * accessible within other local operators.
+ *
+ * @tparam ModelState The model state type
+ */
template<class ModelState>
struct ModelCoefficientMapperBase
{
@@ -58,13 +108,27 @@ struct ModelCoefficientMapperBase
std::map<std::size_t, std::shared_ptr<XView>> _x_view;
std::map<std::size_t, std::shared_ptr<SolutionVector>> _x;
+ /**
+ * @brief Constructs the model coefficient mapper base
+ * @details The operator in the argument refers to the index to which the
+ * coefficient mapper has to use local coefficients (passed as
+ * argument during evaluation) instead of the external
+ * coefficients. The index refers to the index to which the model
+ * states are mapped from.
+ *
+ * @param[in] this_op The operator for local coefficients
+ */
ModelCoefficientMapperBase(std::size_t this_op)
: _this_op(this_op)
{}
+ /**
+ * @brief Updates the given model states.
+ *
+ * @param[in] states The model states
+ */
void update(const std::map<std::size_t, ModelState>& states)
{
-
_x.clear();
_x_view.clear();
_lfs_cache.clear();
@@ -81,6 +145,15 @@ struct ModelCoefficientMapperBase
}
}
+ /**
+ * @brief Bind internal coefficients to an entity
+ * @details This binds the entity to a internal local function space and
+ * local coefficients to the internal coefficients
+ *
+ * @param[in] entity The entity
+ *
+ * @tparam E The entity
+ */
template<class E>
void bind(const E& entity)
{
@@ -95,6 +168,9 @@ struct ModelCoefficientMapperBase
}
}
+ /**
+ * @brief Unbind this coefficients to an entity
+ */
void unbind()
{
for (auto& [op, x_view] : _x_view)
@@ -102,6 +178,14 @@ struct ModelCoefficientMapperBase
}
};
+/**
+ * @brief Constructs the model coefficient mapper
+ * @details This coefficient mapper is designed to one domain local
+ * operators. That is, when the PDELab local function tree only
+ * has depth equal to 1
+ *
+ * @tparam ModelState The model state type
+ */
template<class ModelState>
struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
{
@@ -113,10 +197,21 @@ struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
// component -> (operator,child)
std::unordered_map<std::size_t, std::array<std::size_t, 2>> _comp_mapper;
+ /**
+ * @brief Constructor of the coefficient mapper
+ *
+ * @param[in] this_op The operator for local coefficients
+ */
ModelCoefficientMapper(std::size_t this_op = 0)
: Base(this_op)
{}
+ /**
+ * @brief Constructor of the coefficient mapper
+ *
+ * @param[in] map_operator A vector mapping components to operators
+ * @param[in] this_op The operator for local coefficients
+ */
ModelCoefficientMapper(const std::vector<std::size_t>& map_operator,
std::size_t this_op)
: Base(this_op)
@@ -124,6 +219,12 @@ struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
set_mapper(map_operator);
}
+ /**
+ * @brief Constructor of the coefficient mapper
+ *
+ * @param[in] operator_config The operator configuration
+ * @param[in] this_op The operator for local coefficients
+ */
ModelCoefficientMapper(const ParameterTree& operator_config,
std::size_t this_op)
: Base(this_op)
@@ -141,6 +242,12 @@ struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
set_mapper(map_op);
}
+ /**
+ * @brief Sets a map from component to operator and child in the PDELab
+ * tree
+ *
+ * @param[in] map_op A vector mapping components to operators
+ */
void set_mapper(const std::vector<std::size_t>& map_op)
{
std::size_t max_op = *std::max_element(map_op.begin(), map_op.end());
@@ -151,6 +258,20 @@ struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
std::array<std::size_t, 2>{ map_op[i], comp_child[map_op[i]]++ };
}
+ /**
+ * @brief Evaluates the coefficents for a give component and a degree of
+ * freedom. Either from a local vector x_view_local, or by other
+ * stored coefficients vector.
+ *
+ * @param[in] x_view_local The view on the local coefficients vector
+ * @param[in] comp The component
+ * @param[in] dof The degree of freedom
+ *
+ * @tparam XViewLocal The type fo the local coefficient view
+ *
+ * @return The coefficient value for the requested component and the
+ * degree of freedom
+ */
template<class XViewLocal>
auto operator()(const XViewLocal& x_view_local,
const std::size_t& comp,
@@ -167,6 +288,14 @@ struct ModelCoefficientMapper : public ModelCoefficientMapperBase<ModelState>
}
};
+/**
+ * @brief Constructs the coefficient mapper for multidomain models
+ * @details This coefficient mapper is designed to multidomain domain local
+ * operators. That is, when the PDELab local function tree has
+ * depth equal to 2
+ *
+ * @tparam ModelState The model state type
+ */
template<class ModelState>
struct MultiDomainModelCoefficientMapper
: public ModelCoefficientMapper<ModelState>
@@ -179,6 +308,20 @@ struct MultiDomainModelCoefficientMapper
using Base::Base;
std::size_t _domain;
+ /**
+ * @brief Evaluates the coefficents for a give component and a degree of
+ * freedom. Either from a local vector x_view_local, or by other
+ * stored coefficients vector.
+ *
+ * @param[in] x_view_local The view on the local coefficients vector
+ * @param[in] comp The component
+ * @param[in] dof The degree of freedom
+ *
+ * @tparam XViewLocal The type fo the local coefficient view
+ *
+ * @return The coefficient value for the requested component and the
+ * degree of freedom
+ */
template<class XViewLocal>
auto operator()(const XViewLocal& x_view_local,
const std::size_t& comp,
@@ -195,6 +338,16 @@ struct MultiDomainModelCoefficientMapper
}
}
+ /**
+ * @brief Bind internal coefficients to an entity
+ * @details This procedure sets the first child in the PDELab tree by
+ * finding the domain the entity belongs to, plus further bindings
+ * for the other coefficients
+ *
+ * @param[in] entity The entity
+ *
+ * @tparam E The entity
+ */
template<class E>
void bind(const E& entity)
{
diff --git a/dune/copasi/enum.hh b/dune/copasi/common/enum.hh
similarity index 67%
rename from dune/copasi/enum.hh
rename to dune/copasi/common/enum.hh
index a3035334865d0acf78507875c9cbc9b95ca99fe0..c0e90b30ec8fa8e920ebe62f05e0bef68b85b6f8 100644
--- a/dune/copasi/enum.hh
+++ b/dune/copasi/common/enum.hh
@@ -3,6 +3,9 @@
namespace Dune::Copasi {
+/**
+ * @brief This class describes a model setup policy.
+ */
enum class ModelSetupPolicy
{
None,
@@ -16,11 +19,17 @@ enum class ModelSetupPolicy
All
};
+/**
+ * @brief This class describes an adaptivity policy.
+ */
enum class AdaptivityPolicy
{
None
};
+/**
+ * @brief This class describes a jacobian method.
+ */
enum class JacobianMethod
{
Analytical,
diff --git a/dune/copasi/common/muparser_data_handler.hh b/dune/copasi/common/muparser_data_handler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48955f7f31a324bd13b20f5d3f64c45c6200f895
--- /dev/null
+++ b/dune/copasi/common/muparser_data_handler.hh
@@ -0,0 +1,111 @@
+#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 {
+
+/**
+ * @brief MuParser data handler.
+ * @details MuParser only allows to define functions if they are defined
+ * statically. Hence, this is a trick to make happy muparser about
+ * them. We basically store all possible functions in a static
+ * vector, allow mu parser to call them in a static function. The
+ * trick is that we assig the content of each function at run time.
+ * This is mainly for tiff files, but it can be extended for other
+ * cases.
+ * @todo Add static check to check that T is a callable of two arguments
+ *
+ * @tparam T The type to manage statically. It has to have the
+ * parenthesis operator for two arguments.
+ */
+template<class T>
+struct MuParserDataHandler
+{
+
+ /**
+ * @brief Destroys the object.
+ */
+ ~MuParserDataHandler()
+ {
+ MuParserDataHandler<T>::_names.clear();
+ MuParserDataHandler<T>::_functions.clear();
+ }
+
+ /**
+ * @brief The static functions for MuParser
+ *
+ * @param[in] x The x position in global coordinates
+ * @param[in] y The x position in global coordinates
+ *
+ * @tparam i Index to make each function to be unique
+ *
+ * @return The return value of the wrapped function
+ */
+ template<int i>
+ static double function_wrapper(double x, double y)
+ {
+ assert(MuParserDataHandler<T>::_functions[i]);
+ return (*MuParserDataHandler<T>::_functions[i])(x, y);
+ }
+
+ /**
+ * @brief Adds tiff functions.
+ * @todo Disable this function if T is not a tiff grayscale
+ *
+ * @param[in] data_config The data configuration
+ */
+ 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));
+ }
+ }
+
+ /**
+ * @brief Sets the functions to muParser.
+ *
+ * @param parser The muParser parser
+ *
+ * @tparam max_functions Maximum set of static functions to set to
+ * muParser
+ */
+ template<std::size_t max_functions = 20>
+ void set_functions(mu::Parser& parser)
+ {
+ auto indices =
+ Dune::range(std::integral_constant<std::size_t, max_functions>{});
+ Dune::Hybrid::forEach(indices, [&](auto i) {
+ if (i < _functions.size())
+ parser.DefineFun(_names[i], function_wrapper<i>);
+ });
+ }
+
+ /// vector of strigs descriving the name for each function
+ static std::vector<std::string> _names;
+
+ /// vector of callables to be used for each static function
+ static std::vector<std::shared_ptr<T>> _functions;
+};
+
+/// initialization fo the static vectore of names
+template<class T>
+std::vector<std::string> MuParserDataHandler<T>::_names = {};
+
+/// initialization of vector of callables
+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
diff --git a/dune/copasi/common/parameter_tree_check.hh b/dune/copasi/common/parameter_tree_check.hh
new file mode 100644
index 0000000000000000000000000000000000000000..600023e72d71501d9462e6349ce0e94adaba2f4f
--- /dev/null
+++ b/dune/copasi/common/parameter_tree_check.hh
@@ -0,0 +1,282 @@
+#ifndef DUNE_COPASI_PARAMETER_PARSER_HH
+#define DUNE_COPASI_PARAMETER_PARSER_HH
+
+#include <functional>
+
+#include <dune/common/parametertree.hh>
+
+namespace Dune::Copasi {
+
+bool
+eq(const ParameterTree& config_l, const ParameterTree& config_r)
+{
+ // absolutely not obtimal for big trees!
+ for (auto&& key : config_r.getValueKeys())
+ if (not config_l.hasKey(key))
+ return false;
+
+ for (auto&& key : config_l.getValueKeys())
+ if (not config_r.hasKey(key))
+ return false;
+
+ for (auto&& section : config_r.getSubKeys())
+ if (not eq(config_l.sub(section), config_r.sub(section)))
+ return false;
+
+ for (auto&& section : config_l.getSubKeys())
+ if (not eq(config_l.sub(section), config_r.sub(section)))
+ return false;
+
+ return true;
+}
+
+// merge parameter trees
+// config_l += config_r
+void
+add(ParameterTree& config_l,
+ const ParameterTree& config_r,
+ bool trow_on_override = true)
+{
+ for (auto&& key : config_r.getValueKeys()) {
+ if (config_l.hasKey(key))
+ DUNE_THROW(RangeError,
+ "config file addition failed due to duplucated key: " << key);
+ config_l[key] = config_r[key];
+ }
+ for (auto&& section : config_r.getSubKeys())
+ add(config_l.sub(section), config_r.sub(section));
+}
+
+// diff parameter trees
+// config_l -= config_r
+void
+diff(ParameterTree& config_l, const ParameterTree& config_r)
+{
+ ParameterTree config_l_copy(config_l);
+ config_l = {};
+
+ for (auto&& key : config_l_copy.getValueKeys())
+ if (not config_r.hasKey(key))
+ config_l[key] = config_l_copy[key];
+
+ for (auto&& section : config_l_copy.getSubKeys()) {
+ diff(config_l_copy.sub(section), config_r.sub(section));
+ if (config_l_copy.sub(section).getValueKeys().size() > 0 or
+ config_l_copy.sub(section).getSubKeys().size() >
+ 0) // remove empty sections
+ config_l.sub(section) = config_l_copy.sub(section);
+ }
+}
+
+// get keys in vector
+ParameterTree
+get_keys(const ParameterTree& config, const std::vector<std::string>& keys)
+{
+ ParameterTree new_config;
+ for (auto&& key : keys)
+ new_config[key] = config[key];
+ return new_config;
+}
+
+// get all keys (no section)
+ParameterTree
+get_keys(const ParameterTree& config)
+{
+ return get_keys(config, config.getValueKeys());
+}
+
+// get sections in vector
+ParameterTree
+get_sections(
+ const ParameterTree& config,
+ const std::vector<std::function<ParameterTree(const ParameterTree&)>>&
+ sections)
+{
+ ParameterTree new_config;
+ for (auto&& section : sections)
+ add(new_config, section(config));
+ return new_config;
+}
+
+ParameterTree
+get_grid(const ParameterTree& config)
+{
+ std::vector<std::string> grid_keys{ "file", "initial_level" };
+ auto grid_config = get_keys(config.sub("grid"), grid_keys);
+ // todo : check file is valid
+ // todo : check initial_level is > 0
+
+ ParameterTree new_config;
+ new_config.sub("grid") = grid_config;
+ return new_config;
+}
+
+ParameterTree
+get_initial(const ParameterTree& config)
+{
+ auto initial_config = get_keys(config.sub("initial"));
+ // todo : check that they are math expressions
+ // todo : keys are ordered
+
+ ParameterTree new_config;
+ new_config.sub("initial") = initial_config;
+ return new_config;
+}
+
+// check_var_consistency with respect to initial sections
+ParameterTree
+get_diffusion(const ParameterTree& config, bool check_var_consistency = true)
+{
+ ParameterTree diffusion_config;
+ if (check_var_consistency) {
+ auto initial_config = get_keys(config.sub("initial"));
+ diffusion_config =
+ get_keys(config.sub("diffusion"), initial_config.getValueKeys());
+ } else {
+ diffusion_config = get_keys(config.sub("diffusion"));
+ }
+ // todo : check that they are math expressions
+ // todo : keys are ordered
+
+ ParameterTree new_config;
+ new_config.sub("diffusion") = diffusion_config;
+ return new_config;
+}
+
+// check_var_consistency with respect to initial sections
+ParameterTree
+get_operator(const ParameterTree& config, bool check_var_consistency = true)
+{
+ ParameterTree operator_config;
+ if (check_var_consistency) {
+ auto initial_config = get_keys(config.sub("initial"));
+ operator_config =
+ get_keys(config.sub("operator"), initial_config.getValueKeys());
+ } else {
+ operator_config = get_keys(config.sub("operator"));
+ }
+ // todo : check that they are signed integers
+ // todo : keys are ordered
+
+ ParameterTree new_config;
+ new_config.sub("operator") = operator_config;
+ return new_config;
+}
+
+ParameterTree
+get_jacobian(const ParameterTree& config,
+ std::vector<std::string> base_variables)
+{
+ auto jacobian_config = get_keys(config.sub("jacobian"));
+ std::size_t jac_size = jacobian_config.getValueKeys().size();
+ std::size_t base_size = base_variables.size();
+
+ // it's important that jacobian has the right size and is ordered
+ if (jac_size != base_size * base_size)
+ DUNE_THROW(RangeError, "jacobian section has wrong size");
+
+ std::size_t count(0);
+ for (auto&& var_i : base_variables) {
+ for (auto&& var_j : base_variables) {
+ std::string jac_key = jacobian_config.getValueKeys()[count];
+ auto found_i = jac_key.find(var_i);
+ if (found_i == std::string::npos)
+ DUNE_THROW(RangeError,
+ "Jacobian key '" << jac_key
+ << "' does not contain its base key i '"
+ << var_i << "'");
+ auto found_j = jac_key.find(var_j);
+ if (found_j == std::string::npos)
+ DUNE_THROW(RangeError,
+ "Jacobian key '" << jac_key
+ << "' does not contain its base key j '"
+ << var_j << "'");
+ count++;
+ }
+ }
+
+ // todo : check that they are math expressions
+ // todo : keys are ordered
+
+ ParameterTree new_config;
+ new_config.sub("jacobian") = jacobian_config;
+ return new_config;
+}
+
+ParameterTree
+get_reaction(const ParameterTree& config, bool check_var_consistency = true)
+{
+ ParameterTree reaction_config;
+ if (check_var_consistency) {
+ auto initial_config = get_keys(config.sub("initial"));
+ reaction_config =
+ get_keys(config.sub("reaction"), initial_config.getValueKeys());
+ } else {
+ reaction_config = get_keys(config.sub("reaction"));
+ }
+
+ auto jacobian_config =
+ get_jacobian(config.sub("reaction"), reaction_config.getValueKeys());
+ add(reaction_config, jacobian_config);
+ // todo : check that they are math expressions
+ // todo : keys are ordered
+
+ ParameterTree new_config;
+ new_config.sub("reaction") = reaction_config;
+ return new_config;
+}
+
+ParameterTree
+get_compartment_i(const ParameterTree& config, const std::string& name)
+{
+ std::vector<std::function<ParameterTree(const ParameterTree&)>> sections;
+ sections.push_back(get_initial);
+ sections.push_back([](auto i) { return get_diffusion(i); });
+ sections.push_back([](auto i) { return get_reaction(i); });
+ sections.push_back([](auto i) { return get_operator(i); });
+ auto compart_config = get_sections(config.sub(name), sections);
+
+ ParameterTree new_config;
+ new_config.sub(name) = compart_config;
+ return new_config;
+}
+
+ParameterTree
+get_compartments(const ParameterTree& config)
+{
+ auto compart_config = get_keys(config.sub("compartments"));
+ // todo : check that they are valid keys
+
+ ParameterTree new_config;
+ new_config.sub("compartments") = compart_config;
+ return new_config;
+}
+
+ParameterTree
+get_model(const ParameterTree& config)
+{
+ auto compartments = get_compartments(config.sub("model"));
+
+ std::vector<std::function<ParameterTree(const ParameterTree&)>> sections;
+ for (auto&& compartment : compartments.sub("compartments").getValueKeys())
+ sections.push_back(
+ [&](auto ini) { return get_compartment_i(ini, compartment); });
+
+ auto compartment_i = get_sections(config.sub("model"), sections);
+
+ add(compartments, compartment_i);
+
+ ParameterTree new_config;
+ new_config.sub("model") = compartments;
+
+ // bug in parameter tree. If keys go first, sub section is not well set!
+ std::vector<std::string> model_keys{ "begin_time", "end_time", "time_step" };
+ add(new_config.sub("model"), get_keys(config.sub("model"), model_keys));
+
+ // todo : check that they are valid keys
+ return new_config;
+}
+
+} // namespace Dune::Copasi
+
+#endif // DUNE_COPASI_PARAMETER_PARSER_HH
\ No newline at end of file
diff --git a/dune/copasi/pdelab_expression_adapter.hh b/dune/copasi/common/pdelab_expression_adapter.hh
similarity index 70%
rename from dune/copasi/pdelab_expression_adapter.hh
rename to dune/copasi/common/pdelab_expression_adapter.hh
index 60670113d94b35e456bcaf16b81648ef7ab02c16..1a1b48d3f446ec1df2f42302d7906efb7a440872 100644
--- a/dune/copasi/pdelab_expression_adapter.hh
+++ b/dune/copasi/common/pdelab_expression_adapter.hh
@@ -1,7 +1,7 @@
#ifndef DUNE_COPASI_GRID_FUNCTION_EXPRESSION_ADAPTER_HH
#define DUNE_COPASI_GRID_FUNCTION_EXPRESSION_ADAPTER_HH
-#include <dune/copasi/tiff_grayscale.hh>
+#include <dune/copasi/common/tiff_grayscale.hh>
#include <dune/pdelab/common/function.hh>
@@ -18,6 +18,14 @@
namespace Dune::Copasi {
+/**
+ * @brief Converts an interface to match an expression to a PDELab grid
+ * function.
+ * @details The resulting grid view is only for scalar expressions
+ *
+ * @tparam GV Grid View
+ * @tparam RF Range Field
+ */
template<typename GV, typename RF>
class ExpressionToGridFunctionAdapter
: public PDELab::GridFunctionBase<
@@ -27,8 +35,16 @@ class ExpressionToGridFunctionAdapter
public:
using Traits = PDELab::GridFunctionTraits<GV, RF, 1, FieldVector<RF, 1>>;
- //! construct from grid view
-
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] grid_view The grid view
+ * @param[in] equation The math expression
+ * @param[in] do_compile_parser Bool to compile parser at object
+ * construction
+ * @param[in] other_variables Extra varialbes names to be available in the
+ * expression
+ */
ExpressionToGridFunctionAdapter(const GV& grid_view,
const std::string& equation,
bool do_compile_parser = true,
@@ -67,15 +83,32 @@ public:
_logger.debug("ExpressionToGridFunctionAdapter constructed"_fmt);
}
+ /**
+ * @brief Destroys the object.
+ */
~ExpressionToGridFunctionAdapter()
{
_logger.debug("ExpressionToGridFunctionAdapter deconstructed"_fmt);
}
- //! get a reference to the grid view
+ /**
+ * @brief Gets a reference to the grid view
+ *
+ * @return The grid view.
+ */
inline const GV& getGridView() const { return _gv; }
- //! evaluate extended function on element
+ /**
+ * @brief Evaluates extended function on a element
+ *
+ * @param[in] e Entity to operate with
+ * @param[in] x Local coordinates in the entity
+ * @param y Resulting value
+ *
+ * @tparam E Entity
+ * @tparam D Domain
+ * @tparam R Range
+ */
template<class E, class D, class R>
void evaluate(const E& e, const D& x, R& y) const
{
@@ -91,12 +124,22 @@ public:
}
}
+ /**
+ * @brief Updates the given other values set at construction.
+ * @details The extra values have to have the same order as entred in the
+ * object construction
+ *
+ * @param[in] other_value The other values
+ */
inline void update(const DynamicVector<RF>& other_value)
{
assert(_other_value.size() == other_value.size());
_other_value = other_value;
}
+ /**
+ * @brief Compiles the parser and checks that it is able to be evaluated
+ */
void compile_parser()
{
try {
@@ -116,13 +159,18 @@ public:
return _parser;
}
+ /**
+ * @brief Sets the time.
+ *
+ * @param[in] t The new time
+ */
void set_time(double t) { _time = t; }
private:
- /// Output information on the parser error and throw DUNE exception
/**
- * \param e Exception thrown by the parser
- * \throw IOError (always throws)
+ * @brief Output information on the parser error and throw DUNE exception
+ *
+ * @param[in] e Exception thrown by the parser
*/
void handle_parser_error(const mu::Parser::exception_type& e) const
{
diff --git a/dune/copasi/common/tiff_grayscale.hh b/dune/copasi/common/tiff_grayscale.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8be21255636468880adf596d99d29cbeb7b3ccc0
--- /dev/null
+++ b/dune/copasi/common/tiff_grayscale.hh
@@ -0,0 +1,249 @@
+#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 <memory>
+#include <queue>
+#include <string>
+
+namespace Dune::Copasi {
+
+/**
+ * @brief This class describes a tiff grayscale.
+ * @details This class opens, closes, and quere information from a single
+ * tiff image. Data can be querid using matrix like indices with the
+ * bracket operator or using space coordinates using the parethesis
+ * operator.
+ *
+ * @tparam T Unsigned integer type to represent each pixel information.
+ */
+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");
+
+ /**
+ * @brief This class represents a single row of a tiff grayscale
+ */
+ struct TIFFGrayscaleRow
+ {
+ /**
+ * @brief Constructs a tiff grayscale row
+ *
+ * @param tiff_file The tiff file
+ * @param[in] row The row
+ * @param[in] col_size The column size
+ * @param[in] zero Zero based grayscale?
+ */
+ 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);
+ }
+
+ /**
+ * @brief Array indexer column operator.
+ * @details This operator always scales the gayscale value with its
+ * maximum value. This way, the results is independent of the
+ * template argument T being used
+ *
+ * @param[in] col The column for the current row
+ *
+ * @return The result of the array indexer scaled between 0 and 1.
+ */
+ 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;
+ }
+
+ /**
+ * @brief The size of this row
+ *
+ * @return Number of columns in this row
+ */
+ std::size_t size() const { return static_cast<std::size_t>(_col_size); }
+
+ /**
+ * @brief The current row
+ *
+ * @return Current row
+ */
+ 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:
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] filename The filename
+ * @param[in] max_cache The maximum row cache.
+ */
+ TIFFGrayscale(const std::string& filename, std::size_t max_cache = 8)
+ : _tiff_file(TIFFOpen(filename.c_str(), "r"))
+ , _max_cache(max_cache)
+ {
+ 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);
+ }
+
+ /**
+ * @brief Destroys the object.
+ */
+ ~TIFFGrayscale() { TIFFClose(_tiff_file); }
+
+ /**
+ * @brief Array indexer row operator.
+ * @warning If rows are read concurrently, this object has to be copied
+ *
+ * @param[in] row The row index
+ *
+ * @return The row
+ */
+ const TIFFGrayscaleRow& operator[](T row) const
+ {
+ assert((short)row < _row_size);
+ return cache(row);
+ }
+
+ /**
+ * @brief TIFF value by position call operator.
+ * @details Here we assume that TIFFTAG_RESOLUTIONUNIT has same units as
+ * the grid. Since we never check grid units, we also do not check
+ * tiff units
+ *
+ * @param[in] x The x position in same units as the tiff file
+ * @param[in] y The x position in same units as the tiff file
+ *
+ * @tparam DF Domain Field type
+ *
+ * @return Scaled TIFF value if x and y are in the TIFF domain, 0
+ * otherwise
+ */
+ template<class DF>
+ double operator()(const DF& x, const DF& y)
+ {
+ // 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];
+ }
+
+ /**
+ * @brief The size of rows for this file
+ *
+ * @return Number of rows in this file
+ */
+ std::size_t size() const { return rows(); }
+
+ /**
+ * @brief The size of rows for this file
+ *
+ * @return Number of rows in this file
+ */
+ std::size_t rows() const { return static_cast<std::size_t>(_row_size); }
+
+ /**
+ * @brief The size of cols for this file
+ *
+ * @return Number of cols in this file
+ */
+ std::size_t cols() const { return static_cast<std::size_t>(_col_size); }
+
+private:
+ /**
+ * @brief Cache rows
+ *
+ * @param[in] row The row to be cached
+ *
+ * @return The cached row
+ */
+ 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;
+ const std::size_t _max_cache;
+};
+
+} // namespace Dune::Copasi
+
+#endif // DUNE_COPASI_TIFF_GRAYSCALE_HH
\ No newline at end of file
diff --git a/dune/copasi/finite_element/CMakeLists.txt b/dune/copasi/finite_element/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eac0d20d8ff38c50ca14077b97e6ea1d600c7070
--- /dev/null
+++ b/dune/copasi/finite_element/CMakeLists.txt
@@ -0,0 +1,6 @@
+install(FILES dynamic_local_basis.hh
+ dynamic_local_coefficients.hh
+ dynamic_local_interpolation.hh
+ dynamic_local_finite_element.hh
+ multidomain_local_finite_element_map.hh
+ DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/dune/copasi/finite_element")
\ No newline at end of file
diff --git a/dune/copasi/dynamic_local_basis.hh b/dune/copasi/finite_element/dynamic_local_basis.hh
similarity index 50%
rename from dune/copasi/dynamic_local_basis.hh
rename to dune/copasi/finite_element/dynamic_local_basis.hh
index 169e672481f8565198ead7c81ab6307cce07dfef..abcc7b7afeecaf3bd3fa2a9e5b1c5db8f80ea484 100644
--- a/dune/copasi/dynamic_local_basis.hh
+++ b/dune/copasi/finite_element/dynamic_local_basis.hh
@@ -9,12 +9,23 @@
namespace Dune::Copasi {
+/**
+ * @brief This class describes dynamic power local basis.
+ *
+ * @tparam Basis The base local basis
+ */
template<class Basis>
class DynamicPowerLocalBasis
{
public:
using Traits = typename Basis::Traits;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] basis The base local basis
+ * @param[in] power_size The power size
+ */
DynamicPowerLocalBasis(const Basis& basis, std::size_t power_size)
: _power_size(power_size)
, _basis(basis)
@@ -23,22 +34,53 @@ public:
assert(_power_size >= 0);
}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] basis The base local basis
+ */
DynamicPowerLocalBasis(const Basis& basis)
: DynamicPowerLocalBasis(basis, 1)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] power_size The power size
+ *
+ * @tparam <unnamed> Internal use to allow default constructible base
+ * local basis
+ */
template<class = std::enable_if_t<std::is_default_constructible_v<Basis>>>
DynamicPowerLocalBasis(std::size_t power_size)
: DynamicPowerLocalBasis(Basis{}, power_size)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @tparam <unnamed> Internal use to allow default constructible base
+ * local basis
+ */
template<class = std::enable_if_t<std::is_default_constructible_v<Basis>>>
DynamicPowerLocalBasis()
: DynamicPowerLocalBasis(1)
{}
+ /**
+ * @brief Returns the number of shape functions for this local finite
+ * element
+ *
+ * @return The number of shape functions
+ */
unsigned int size() const { return _power_size * _basis.size(); }
+ /**
+ * @brief Evaluation function for a given local coordinate
+ *
+ * @param[in] in The local coordinates
+ * @param out The value of each shape function
+ */
inline void evaluateFunction(
const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const
@@ -47,6 +89,12 @@ public:
populate_output(f, in, out);
}
+ /**
+ * @brief Jacobian evaluation function for a given local coordinate
+ *
+ * @param[in] in The local coordinates
+ * @param out The value of the jacobian for each shape function
+ */
inline void evaluateJacobian(
const typename Traits::DomainType& in,
std::vector<typename Traits::JacobianType>& out) const
@@ -55,6 +103,18 @@ public:
populate_output(f, in, out);
}
+ /**
+ * @brief The partial derivate of the shape functions
+ * @warning This function is here only fulfill the interface requirements.
+ * This was never used nor tested.
+ *
+ * @param[in] order The order
+ * @param[in] in The local coordinates
+ * @param out The value of the partial derivate for each shape
+ * function
+ *
+ * @tparam dim The grid entity dimension
+ */
template<int dim>
void partial(const std::array<unsigned int, dim>& order,
const typename Traits::DomainType& in,
@@ -64,9 +124,27 @@ public:
populate_output(f, in, out);
}
+ /**
+ * @brief Polynomal order of the shape functions
+ *
+ * @return The polynomal order of the shape functions
+ */
unsigned int order() const { return _basis.order(); }
private:
+ /**
+ * @brief Helper class to polulate results for interface functions
+ *
+ * @param[in] f A function that evaluates a interface function for an
+ * individual local basis function
+ * @param[in] in The local coordinates
+ * @param out The result of the evaluation of f over the power local
+ * basis function
+ *
+ * @tparam F The function type
+ * @tparam In The local coordinates type
+ * @tparam Out The result populated type
+ */
template<class F, class In, class Out>
void populate_output(const F& f, const In& in, Out& out) const
{
diff --git a/dune/copasi/dynamic_local_coefficients.hh b/dune/copasi/finite_element/dynamic_local_coefficients.hh
similarity index 68%
rename from dune/copasi/dynamic_local_coefficients.hh
rename to dune/copasi/finite_element/dynamic_local_coefficients.hh
index 47ac7d2a5ccbcbb5377fe3b106243aedd695ec02..b1fb8a701bd47cbb347e48efd331eb7ad674fa09 100644
--- a/dune/copasi/dynamic_local_coefficients.hh
+++ b/dune/copasi/finite_element/dynamic_local_coefficients.hh
@@ -11,10 +11,21 @@
namespace Dune::Copasi {
+/**
+ * @brief This class describes dynamic power local coefficients.
+ *
+ * @tparam Coefficients The base local coefficients
+ */
template<class Coefficients>
class DynamicPowerLocalCoefficients
{
public:
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] coefficients The base coefficients
+ * @param[in] power_size The power size
+ */
DynamicPowerLocalCoefficients(const Coefficients& coefficients,
std::size_t power_size)
: _size(power_size * coefficients.size())
@@ -53,24 +64,55 @@ public:
assert(_local_keys.size() == size());
}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] coefficients The base coefficients
+ */
DynamicPowerLocalCoefficients(const Coefficients& coefficients)
: DynamicPowerLocalCoefficients(coefficients, 1)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] power_size The power size
+ *
+ * @tparam <unnamed> Internal use to allow default constructible base
+ * local coefficients
+ */
template<
class = std::enable_if_t<std::is_default_constructible_v<Coefficients>>>
DynamicPowerLocalCoefficients(std::size_t power_size)
: DynamicPowerLocalCoefficients(Coefficients{}, power_size)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @tparam <unnamed> Internal use to allow default constructible base
+ * local coefficients
+ */
template<
class = std::enable_if_t<std::is_default_constructible_v<Coefficients>>>
DynamicPowerLocalCoefficients()
: DynamicPowerLocalCoefficients(1)
{}
+ /**
+ * @brief Returns the number of coefficients
+ *
+ * @return The size
+ */
std::size_t size() const { return _size; }
+ /**
+ * @brief Returns a local key for each degree of freedom
+ *
+ * @param[in] i The degree of freedom
+ *
+ * @return The local key for the given degree of freedom
+ */
const LocalKey& localKey(std::size_t i) const { return _local_keys[i]; }
private:
diff --git a/dune/copasi/dynamic_local_finite_element.hh b/dune/copasi/finite_element/dynamic_local_finite_element.hh
similarity index 59%
rename from dune/copasi/dynamic_local_finite_element.hh
rename to dune/copasi/finite_element/dynamic_local_finite_element.hh
index ccc8a78075bd67464347efab0bf5829dce532b35..fc92602eb276f183e1e5ec39a810cdb47f897d91 100644
--- a/dune/copasi/dynamic_local_finite_element.hh
+++ b/dune/copasi/finite_element/dynamic_local_finite_element.hh
@@ -1,9 +1,9 @@
#ifndef DUNE_COPASI_DYNAMIC_LOCAL_FINITE_ELEMENT_HH
#define DUNE_COPASI_DYNAMIC_LOCAL_FINITE_ELEMENT_HH
-#include <dune/copasi/dynamic_local_basis.hh>
-#include <dune/copasi/dynamic_local_coefficients.hh>
-#include <dune/copasi/dynamic_local_interpolation.hh>
+#include <dune/copasi/finite_element/dynamic_local_basis.hh>
+#include <dune/copasi/finite_element/dynamic_local_coefficients.hh>
+#include <dune/copasi/finite_element/dynamic_local_interpolation.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
@@ -11,6 +11,11 @@
namespace Dune::Copasi {
+/**
+ * @brief This class describes a dynamic power local finite element.
+ *
+ * @tparam LocalFiniteElement The base local finite element
+ */
template<class LocalFiniteElement>
class DynamicPowerLocalFiniteElement
{
@@ -30,6 +35,14 @@ public:
using Traits =
LocalFiniteElementTraits<LocalBasis, LocalCoefficients, LocalInterpolation>;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] finite_element The base local finite element
+ * @param[in] local_basis The base local basis
+ * @param[in] local_coeficcients The base local coeficcients
+ * @param[in] local_interpolation The base local interpolation
+ */
DynamicPowerLocalFiniteElement(const LocalFiniteElement& finite_element,
const LocalBasis& local_basis,
const LocalCoefficients& local_coeficcients,
@@ -40,6 +53,15 @@ public:
, _interpolation(local_interpolation)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] finite_element The base local finite element
+ * @param[in] power_size The power size
+ * @param[in] local_basis The base local basis
+ * @param[in] local_coeficcients The base local coeficcients
+ * @param[in] local_interpolation The base local interpolation
+ */
DynamicPowerLocalFiniteElement(const LocalFiniteElement& finite_element,
std::size_t power_size)
: _finite_element(finite_element)
@@ -48,24 +70,54 @@ public:
, _interpolation(_finite_element.localInterpolation(), power_size)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] power_size The power size
+ */
DynamicPowerLocalFiniteElement(std::size_t power_size)
: DynamicPowerLocalFiniteElement(LocalFiniteElement{}, power_size)
{}
+ /**
+ * @brief Returns the power local basis
+ *
+ * @return The power local basis
+ */
const typename Traits::LocalBasisType& localBasis() const { return _basis; }
+ /**
+ * @brief Returns the power local coefficents
+ *
+ * @return The power local coefficients
+ */
const typename Traits::LocalCoefficientsType& localCoefficients() const
{
return _coefficients;
}
+ /**
+ * @brief Returns the power local interpolation
+ *
+ * @return The power local interpolation
+ */
const typename Traits::LocalInterpolationType& localInterpolation() const
{
return _interpolation;
}
+ /**
+ * @brief Returns the number of degrees of freedom for this element
+ *
+ * @return The size
+ */
unsigned int size() const { return _coefficients.size(); }
+ /**
+ * @brief Returns the geometry type for this finite element
+ *
+ * @return The geometry type
+ */
GeometryType type() const { return _finite_element.type(); }
private:
diff --git a/dune/copasi/dynamic_local_interpolation.hh b/dune/copasi/finite_element/dynamic_local_interpolation.hh
similarity index 69%
rename from dune/copasi/dynamic_local_interpolation.hh
rename to dune/copasi/finite_element/dynamic_local_interpolation.hh
index f6a0c5f47f043d6bd06bfae580c9daa99cd8f7e2..cfe518745494aeaae9101c22d47a0611b13c33d1 100644
--- a/dune/copasi/dynamic_local_interpolation.hh
+++ b/dune/copasi/finite_element/dynamic_local_interpolation.hh
@@ -6,10 +6,16 @@
#include <dune/common/classname.hh>
#include <dune/common/dynvector.hh>
#include <dune/common/fvector.hh>
+
#include <iostream>
namespace Dune::Copasi {
+/**
+ * @brief This class describes a dynamic power local interpolation.
+ *
+ * @tparam Interpolation The base interpolation
+ */
template<class Interpolation>
class DynamicPowerLocalInterpolation
{
@@ -22,22 +28,52 @@ public:
assert(_power_size >= 0);
}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] interpolation The base interpolation
+ */
DynamicPowerLocalInterpolation(const Interpolation& interpolation)
: DynamicPowerLocalInterpolation(interpolation, 1)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] power_size The power size
+ *
+ * @tparam <unnamed> Internal template argument to allow default
+ * constructible base interpolations
+ */
template<
class = std::enable_if_t<std::is_default_constructible_v<Interpolation>>>
DynamicPowerLocalInterpolation(std::size_t power_size)
: DynamicPowerLocalInterpolation(Interpolation{}, power_size)
{}
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @tparam <unnamed> Internal template argument to allow default
+ * constructible base interpolations
+ */
template<
class = std::enable_if_t<std::is_default_constructible_v<Interpolation>>>
DynamicPowerLocalInterpolation()
: DynamicPowerLocalInterpolation(1)
{}
+ /**
+ * @brief Interpolation method
+ *
+ * @param[in] f Function to be interpolated. If its return type is a
+ * dynamic vector, it has to have a size equal or greater
+ * than the power size set to this ibject
+ * @param out The vector of coefficients for the local finite element
+ *
+ * @tparam F The function type
+ * @tparam C The coefficents type
+ */
template<typename F, typename C>
void interpolate(const F& f, std::vector<C>& out) const
{
diff --git a/dune/copasi/multidomain_local_finite_element_map.hh b/dune/copasi/finite_element/multidomain_local_finite_element_map.hh
similarity index 52%
rename from dune/copasi/multidomain_local_finite_element_map.hh
rename to dune/copasi/finite_element/multidomain_local_finite_element_map.hh
index 1d62a11c86cd3eb5508db45bab9ce381a5e08838..ece64b44b2c03fd5e3d9b6d0f3ee7a534fd97827 100644
--- a/dune/copasi/multidomain_local_finite_element_map.hh
+++ b/dune/copasi/finite_element/multidomain_local_finite_element_map.hh
@@ -1,12 +1,26 @@
#ifndef DUNE_COPASI_MULTIDOMAIN_LOCAL_FINITE_ELEMENT_MAP_HH
#define DUNE_COPASI_MULTIDOMAIN_LOCAL_FINITE_ELEMENT_MAP_HH
-#include <dune/pdelab/finiteelementmap/finiteelementmap.hh>
+#include <dune/copasi/finite_element/dynamic_local_finite_element.hh>
-#include <dune/copasi/dynamic_local_finite_element.hh>
+#include <dune/pdelab/finiteelementmap/finiteelementmap.hh>
namespace Dune::Copasi {
+/**
+ * @brief This class describes a multi domain local finite element map.
+ * @details This class wrapps a usual PDELab finite element map into a
+ * dynamic power finite element map. If the entity to be map does
+ * not belong to the grid view, it will return a 0 power element map
+ * that turns to return a finite element with no degrees of freedom.
+ * This behaviour is useful when the pdelab machinary is operating
+ * in the whole grid (e.g. multidomain grid) but you want to have
+ * different finite elements per sub domain.
+ *
+ * @tparam FiniteElementMap The original finite element map to wrap
+ * @tparam GridView The grid view where the finite element will be
+ * not zero
+ */
template<class FiniteElementMap, class GridView>
class MultiDomainLocalFiniteElementMap
: public PDELab::LocalFiniteElementMapInterface<
@@ -20,6 +34,17 @@ class MultiDomainLocalFiniteElementMap
public:
static const int dimension = FiniteElementMap::dimension;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] grid_view The grid view (may be a sub domain)
+ * @param[in] fem The finite element map to wrap
+ * @param[in] base_finite_element The base finite element of the finite
+ * element map
+ * @param[in] power_size The power size for the sub domain part
+ * (this is always 0 outside of the
+ * subdomain)
+ */
MultiDomainLocalFiniteElementMap(GridView grid_view,
FiniteElementMap fem,
BaseFiniteElement base_finite_element,
@@ -31,6 +56,19 @@ public:
, _fe_null(base_finite_element, 0)
{}
+ /**
+ * @brief Constructs a new instance.
+ * @details This constructor is only available if the base finite element
+ * is default constructiible
+ *
+ * @param[in] grid_view The grid view
+ * @param[in] fem The finite element map to wrap
+ * @param[in] power_size The power size for the sub domain part (this is
+ * always 0 outside of the subdomain)
+ *
+ * @tparam <unnamed> Template helper to allow default construction of
+ * base finite element constructor. Internal use only
+ */
template<
typename =
std::enable_if_t<std::is_default_constructible_v<BaseFiniteElement>, int>>
@@ -43,8 +81,22 @@ public:
power_size)
{}
+ /**
+ * @brief Destroys the object.
+ */
~MultiDomainLocalFiniteElementMap() { delete _fe_cache; }
+ /**
+ * @brief Searches for the finite element for entity e.
+ * @todo Cache more than one finite element
+ *
+ * @param[in] e The entity
+ *
+ * @tparam EntityType The entity
+ *
+ * @return A dynamic power local finite element
+ * @DynamicPowerLocalFiniteElement
+ */
template<class EntityType>
const FiniteElement& find(const EntityType& e) const
{
@@ -66,12 +118,15 @@ public:
.contains(sub_domain);
if (not in_grid_view)
return _fe_null;
+ } else {
+ DUNE_THROW(NotImplemented,
+ "Method not implemented for subdomain entites");
}
}
auto base_fe = _fem.find(e);
+
// cache the last used base finite element
- if (_base_fe_cache != &base_fe) // TODO: Cache more than one fe
- {
+ if (_base_fe_cache != &base_fe) {
_base_fe_cache = &base_fe;
if (_fe_cache != NULL)
delete _fe_cache;
@@ -80,18 +135,42 @@ public:
return *_fe_cache;
}
+ /**
+ * @brief Returns true if this finite element map has a fixed size
+ *
+ * @return Always false for this type
+ */
bool fixedSize() const { return false; }
+ /**
+ * @brief Size for a given geometry type
+ *
+ * @param[in] gt The geometry type
+ *
+ * @return The size for the given geometry type
+ */
std::size_t size(GeometryType gt) const
{
return _power_size * _fem.size(gt);
}
+ /**
+ * @brief Describes wheter a codim has associated degrees of freedom
+ *
+ * @param[in] codim The codim
+ *
+ * @return True if codim has dregees of freedom
+ */
bool hasDOFs(int codim) const
{
return (_power_size != 0) && _fem.hasDOFs(codim);
}
+ /**
+ * @brief Maximum local size for all finite elements
+ *
+ * @return The maximim size
+ */
std::size_t maxLocalSize() const { return _power_size * _fem.maxLocalSize(); }
private:
diff --git a/dune/copasi/grid/CMakeLists.txt b/dune/copasi/grid/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4c65b42015765434df6bfa6a4294f1033f45d438
--- /dev/null
+++ b/dune/copasi/grid/CMakeLists.txt
@@ -0,0 +1,3 @@
+install(FILES multidomain_gmsh_reader.hh
+ multidomain_entity_transformation.hh
+ DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/dune/copasi/grid")
\ No newline at end of file
diff --git a/dune/copasi/multidomain_entity_transformation.hh b/dune/copasi/grid/multidomain_entity_transformation.hh
similarity index 57%
rename from dune/copasi/multidomain_entity_transformation.hh
rename to dune/copasi/grid/multidomain_entity_transformation.hh
index ef0bee2b3780d605a537040c967ec4432681f2b0..6e1077462597f34fc12999be1998e2d8613b574b 100644
--- a/dune/copasi/multidomain_entity_transformation.hh
+++ b/dune/copasi/grid/multidomain_entity_transformation.hh
@@ -5,6 +5,17 @@
namespace Dune::Copasi {
+/**
+ * @brief Multidomain entity transformation
+ * @details This class transform an entity of a multidomain grid to a
+ * subdomain grid entity. Equivalent to the lambda function:
+ * @code{.cpp}
+ * auto entity_transform = [&](auto e){return
+ * _grid->multiDomainEntity(entity);};
+ * @endcode
+ *
+ * @tparam Grid The grid
+ */
template<class Grid>
struct MultiDomainEntityTransformation
{
@@ -14,15 +25,28 @@ struct MultiDomainEntityTransformation
using SubDomainEntity =
typename Grid::SubDomainGrid::Traits::template Codim<0>::Entity;
+ /**
+ * @brief Constructor.
+ *
+ * @param[in] grid A pointer to the grid
+ */
MultiDomainEntityTransformation(std::shared_ptr<const Grid> grid)
: _grid(grid)
{}
+ /**
+ * @brief Function call operator.
+ *
+ * @param[in] entity The subdomain entity
+ *
+ * @return The multidomain entity
+ */
const MultiDomainEntity& operator()(const SubDomainEntity& entity)
{
return _grid->multiDomainEntity(entity);
}
+private:
const std::shared_ptr<const Grid> _grid;
};
diff --git a/dune/copasi/gmsh_reader.hh b/dune/copasi/grid/multidomain_gmsh_reader.hh
similarity index 62%
rename from dune/copasi/gmsh_reader.hh
rename to dune/copasi/grid/multidomain_gmsh_reader.hh
index c5996759bbf8cbf2822a7d46fe5034b5e84770fb..2a17b1701d549aa0e4ef8ddb692e26e292b71072 100644
--- a/dune/copasi/gmsh_reader.hh
+++ b/dune/copasi/grid/multidomain_gmsh_reader.hh
@@ -1,26 +1,63 @@
#ifndef DUNE_COPASI_GMSH_READER_HH
#define DUNE_COPASI_GMSH_READER_HH
+#include <dune/copasi/concepts/grid.hh>
+
#include <dune/grid/common/gridinfo.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/multidomaingrid.hh>
#include <dune/logging/logging.hh>
+#include <dune/common/typetraits.hh>
+
#include <algorithm>
#include <type_traits>
namespace Dune::Copasi {
+/**
+ * @brief This class describes a multi domain gmsh reader.
+ *
+ * @tparam Grid The grid
+ */
template<class Grid>
-class GmshReader;
+class MultiDomainGmshReader
+{
+ static_assert(Dune::Copasi::Concept::isMultiDomainGrid<Grid>(),
+ "Grid is not multidomain grid");
+
+ // this is due the fact that multidomain grids do not export its traits type.
+ // And, since it is needed for the some grids construction, the only way is by
+ // template specialization.
+ static_assert(AlwaysFalse<Grid>::value,
+ "Not implemented: This class only accepts multidomain grids "
+ "from the dune-multidomain grid module");
+};
+/**
+ * @brief This class describes a multi domain gmsh reader.
+ *
+ * @tparam HostGrid The host grid of the multidomain grid
+ * @tparam MDGTraits Traits of the multidomain grid (this is different to
+ * the grid interface traits!)
+ */
template<class HostGrid, class MDGTraits>
-class GmshReader<Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>>
+class MultiDomainGmshReader<Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>>
{
public:
using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+ /**
+ * @brief Reads a grid out of a file name and a configuration file
+ *
+ * @param[in] fileName The gmsh file name
+ * @param[in] config The configuration file
+ * @param[in] insertBoundarySegments Bool to include boundary segments while
+ * reading file
+ *
+ * @return A pair containing pointers to the host and the multidomain grid
+ */
static auto read(const std::string& fileName,
ParameterTree config,
bool insertBoundarySegments = true)
diff --git a/dune/copasi/grid_function_writer.hh b/dune/copasi/grid_function_writer.hh
deleted file mode 100644
index d155566e9d7ec922555a5eeb87b5bac1a4a1ae74..0000000000000000000000000000000000000000
--- a/dune/copasi/grid_function_writer.hh
+++ /dev/null
@@ -1,240 +0,0 @@
-#ifndef DUNE_COPASI_GRID_FUNCTION_WRITER_HH
-#define DUNE_COPASI_GRID_FUNCTION_WRITER_HH
-
-#include <dune/copasi/concepts/grid.hh>
-#include <dune/copasi/concepts/typetree.hh>
-
-#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
-
-#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
-#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/parametertree.hh>
-
-#include <memory>
-#include <string>
-#include <type_traits>
-
-namespace Dune::Copasi {
-
-template<typename Data,
- typename GV = typename Data::LFS::Traits::GridView,
- typename RF =
- typename BasisInterfaceSwitch<typename FiniteElementInterfaceSwitch<
- typename Data::LFS::Traits::FiniteElement>::Basis>::RangeField,
- int dimRange =
- BasisInterfaceSwitch<typename FiniteElementInterfaceSwitch<
- typename Data::LFS::Traits::FiniteElement>::Basis>::dimRange,
- typename Range =
- typename BasisInterfaceSwitch<typename FiniteElementInterfaceSwitch<
- typename Data::LFS::Traits::FiniteElement>::Basis>::Range>
-class DiscreteGridFunction
- : public PDELab::GridFunctionBase<
- PDELab::GridFunctionTraits<GV, RF, dimRange, Range>,
- DiscreteGridFunction<Data, GV, RF, dimRange, Range>>
-{
- typedef PDELab::GridFunctionBase<
- PDELab::GridFunctionTraits<GV, RF, dimRange, Range>,
- DiscreteGridFunction<Data, GV, RF, dimRange, Range>>
- BaseT;
-
-public:
- typedef typename BaseT::Traits Traits;
-
- DiscreteGridFunction(const shared_ptr<Data>& data,
- GV grid_view,
- std::size_t id_begin,
- std::size_t id_end)
- : _grid_view(grid_view)
- , _data(data)
- , _basis(data->_lfs.maxSize())
- , _id_begin(id_begin)
- , _id_end(id_end)
- {}
-
- // Evaluate
- void evaluate(const typename Traits::ElementType& e,
- const typename Traits::DomainType& x,
- typename Traits::RangeType& y) const
- {
- y = 0;
-
- using LFS = decltype(_data->_lfs);
-
- if constexpr (Concept::isSubDomainGrid<typename GV::Grid>()) {
- _data->bind(_grid_view.grid().multiDomainGrid().multiDomainEntity(e));
-
- if constexpr (Concept::isTypeTreeNode<LFS>()) {
- int subdomain = _grid_view.grid().domain();
- auto lfs_child = _data->_lfs.child(subdomain);
- typedef FiniteElementInterfaceSwitch<typename decltype(
- lfs_child)::Traits::FiniteElement>
- FESwitch;
-
- FESwitch::basis(lfs_child.finiteElement()).evaluateFunction(x, _basis);
- for (unsigned int i = _id_begin; i < _id_end; i++)
- y.axpy(_data->_x_local(lfs_child, i), _basis[i]);
- } else {
- static_assert(AlwaysTrue<GV>::value);
- }
- } else {
- _data->bind(e);
-
- typedef FiniteElementInterfaceSwitch<typename LFS::Traits::FiniteElement>
- FESwitch;
-
- FESwitch::basis(_data->_lfs.finiteElement()).evaluateFunction(x, _basis);
-
- for (unsigned int i = _id_begin; i < _id_end; i++)
- y.axpy(_data->_x_local(_data->_lfs, i), _basis[i]);
- }
- }
-
- //! get a reference to the GridView
- const typename Traits::GridViewType& gridView() const { return _grid_view; }
-
-private:
- GV _grid_view;
- const shared_ptr<Data> _data;
- mutable std::vector<typename Traits::RangeType> _basis;
- const std::size_t _id_begin, _id_end;
-};
-
-/**
- * @brief Converts PDELab grid function to match a dune-function.
- * @note This adapter differs from the one provided in PDELab in that this
- * one has the current interface of dune-grid to write vtk files.
- * The main difference using this adapter is that it allows to write
- * more than 3 components.
- *
- * @tparam GF PDELab grid function type.
- */
-template<class GF>
-class VTKGridFunctionAdapter
-{
- using GridFunction = typename std::decay<GF>::type;
- using GridView = typename GF::Traits::GridViewType;
- using Entity = typename GridView::template Codim<0>::Entity;
- using Range = typename GF::Traits::RangeType;
- using Domain = typename GF::Traits::DomainType;
-
-public:
- /**
- * @brief Constructs the object
- *
- * @param[in] gf pointer to a constant grid function
- */
- VTKGridFunctionAdapter(std::shared_ptr<const GF> gf)
- : _entity(nullptr)
- , _gf(gf)
- {}
-
- /**
- * @brief Binds an entity to the grid function
- *
- * @param[in] e Entity
- */
- void bind(const Entity& e) { _entity = &e; }
-
- /**
- * @brief Unbinds the previously binded entity
- */
- void unbind() { _entity = nullptr; }
-
- /**
- * @brief Evaluates the grid function at a given coordinate
- *
- * @param[in] x Local coordinate with respect to the last binded entity
- */
- Range operator()(const Domain& x) const
- {
- Range y;
- assert(_entity);
- _gf->evaluate(*_entity, x, y);
- return y;
- }
-
-private:
- const Entity* _entity;
- std::shared_ptr<const GF> _gf;
-};
-
-/**
- * @brief Class for grid fucntion vtk sequence writer. This class works
- * exacly the same way as a VTKSequence writer but it receives
- * directly smart pointers to grid functions
- *
- * @tparam GV GridView of the dune grid
- */
-template<class GV>
-class GridFunctionVTKSequenceWriter : public Dune::VTKSequenceWriter<GV>
-{
-public:
- // export constructors of vtk sequence writer
- using Dune::VTKSequenceWriter<GV>::VTKSequenceWriter;
-
- ~GridFunctionVTKSequenceWriter() { this->clear(); }
- /**
- * @brief Adds a vertex data.
- *
- * @param[in] gf_ptr A pointer to a grid function of the type GF
- * @param[in] name Name of the variable
- *
- * @tparam GF Type of the grid function
- */
- template<class GF>
- void addVertexData(std::shared_ptr<const GF> gf_ptr, std::string name)
- {
- static_assert(std::is_same<typename GF::Traits::GridViewType, GV>::value,
- "GridFunctionVTKSequenceWriter only works with only one "
- "GridView type.");
-
- VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
-
- const int vtk_ncomps = GF::Traits::dimRange;
- const int dim = GF::Traits::dimDomain;
- auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
- : VTK::FieldInfo::Type::scalar;
- auto vtk_info = VTK::FieldInfo(name, vtk_type, vtk_ncomps);
-
- this->vtkWriter()->addVertexData(vtk_gf, vtk_info);
- }
-
- /**
- * @brief Adds a cell data.
- *
- * @param[in] gf_ptr A pointer to a grid function of the type GF
- * @param[in] name Name of the variable
- *
- * @tparam GF Type of the grid function
- */
- template<class GF>
- void addCellData(std::shared_ptr<const GF> gf_ptr, std::string name)
- {
- static_assert(std::is_same<typename GF::Traits::GridViewType, GV>::value,
- "GridFunctionVTKSequenceWriter only works with only one "
- "GridView type.");
-
- VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
-
- const int vtk_ncomps = GF::Traits::dimRange;
- const int dim = GF::Traits::dimDomain;
- auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
- : VTK::FieldInfo::Type::scalar;
- auto vtk_info = VTK::FieldInfo(name, vtk_type, vtk_ncomps);
-
- this->vtkWriter()->addCellData(vtk_gf, vtk_info);
- }
-
- /**
- * @brief Clear all the attached VTKFunctions.
- * @details This is necessary for each iteration if VTKFunctions are
- * managed internaly with pointers instead of references
- */
- void clear() { this->vtkWriter()->clear(); }
-};
-
-} // namespace Dune::Copasi
-
-#endif // DUNE_COPASI_GRID_FUNCTION_WRITER_HH
\ No newline at end of file
diff --git a/dune/copasi/model/CMakeLists.txt b/dune/copasi/model/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7d137ca6a2df6a58cf4325ae725868d0d6a5c89f
--- /dev/null
+++ b/dune/copasi/model/CMakeLists.txt
@@ -0,0 +1,9 @@
+install(FILES base.hh
+ multidomain_local_operator.hh
+ local_operator.hh
+ diffusion_reaction.cc
+ diffusion_reaction.hh
+ multidomain_diffusion_reaction.cc
+ multidomain_diffusion_reaction.hh
+ state.hh
+ DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/dune/copasi/model")
\ No newline at end of file
diff --git a/dune/copasi/model_base.cc b/dune/copasi/model/base.cc
similarity index 98%
rename from dune/copasi/model_base.cc
rename to dune/copasi/model/base.cc
index 8f2410a4bdb69dae0d70c41ed1e0486f720ddea6..bec0aafe09877343d3d01ba19419f752bd35f1f4 100644
--- a/dune/copasi/model_base.cc
+++ b/dune/copasi/model/base.cc
@@ -2,7 +2,7 @@
#include "config.h"
#endif
-#include <dune/copasi/model_base.hh>
+#include <dune/copasi/model/base.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
diff --git a/dune/copasi/model_base.hh b/dune/copasi/model/base.hh
similarity index 98%
rename from dune/copasi/model_base.hh
rename to dune/copasi/model/base.hh
index fa45f54732526124ef762c9aaa83d36b99a8b83c..a0ff199970e5e9dee7ab6f622ef59a844cbaaa98 100644
--- a/dune/copasi/model_base.hh
+++ b/dune/copasi/model/base.hh
@@ -1,7 +1,7 @@
#ifndef DUNE_COPASI_MODEL_BASE_HH
#define DUNE_COPASI_MODEL_BASE_HH
-#include <dune/copasi/enum.hh>
+#include <dune/copasi/common/enum.hh>
#include <dune/logging/logging.hh>
diff --git a/dune/copasi/model_diffusion_reaction.cc b/dune/copasi/model/diffusion_reaction.cc
similarity index 89%
rename from dune/copasi/model_diffusion_reaction.cc
rename to dune/copasi/model/diffusion_reaction.cc
index a8f0ef4a1c3d69c22798a8d5860bc9689de6825c..770e4464b232567441e46314db9f8a1aca221a50 100644
--- a/dune/copasi/model_diffusion_reaction.cc
+++ b/dune/copasi/model/diffusion_reaction.cc
@@ -7,9 +7,8 @@
* file but a header which has to be included when compiling.
*/
-#include <dune/copasi/model_diffusion_reaction.hh>
-#include <dune/copasi/pdelab_callable_adapter.hh>
-#include <dune/copasi/pdelab_expression_adapter.hh>
+#include <dune/copasi/common/pdelab_expression_adapter.hh>
+#include <dune/copasi/model/diffusion_reaction.hh>
#include <dune/pdelab/function/callableadapter.hh>
#include <dune/pdelab/gridfunctionspace/vtk.hh>
@@ -21,11 +20,11 @@
namespace Dune::Copasi {
template<class Traits>
-ModelDiffusionReaction<Traits>::
- ModelDiffusionReaction(std::shared_ptr<Grid> grid,
- const Dune::ParameterTree& config,
- GV grid_view,
- ModelSetupPolicy setup_policy)
+ModelDiffusionReaction<Traits>::ModelDiffusionReaction(
+ std::shared_ptr<Grid> grid,
+ const Dune::ParameterTree& config,
+ GV grid_view,
+ ModelSetupPolicy setup_policy)
: ModelBase(config)
, _solver_logger(Logging::Logging::componentLogger(config, "solver"))
, _components(config.sub("reaction").getValueKeys().size())
@@ -48,8 +47,7 @@ ModelDiffusionReaction<Traits>::
}
template<class Traits>
-ModelDiffusionReaction<Traits>::
- ~ModelDiffusionReaction()
+ModelDiffusionReaction<Traits>::~ModelDiffusionReaction()
{
_logger.debug("ModelDiffusionReaction deconstructed"_fmt);
}
@@ -87,8 +85,7 @@ ModelDiffusionReaction<Traits>::step()
template<class Traits>
void
-ModelDiffusionReaction<Traits>::suggest_timestep(
- double dt)
+ModelDiffusionReaction<Traits>::suggest_timestep(double dt)
{
DUNE_THROW(NotImplemented, "");
}
@@ -131,8 +128,8 @@ ModelDiffusionReaction<Traits>::suggest_timestep(
template<class Traits>
auto
-ModelDiffusionReaction<Traits>::
- setup_component_grid_function_space(std::string name) const
+ModelDiffusionReaction<Traits>::setup_component_grid_function_space(
+ std::string name) const
{
_logger.trace("create a finite element map"_fmt);
BaseFEM base_fem(_grid->leafGridView());
@@ -140,8 +137,7 @@ ModelDiffusionReaction<Traits>::
_logger.trace("setup grid function space for component {}"_fmt, name);
const ES entity_set(_grid->leafGridView());
- auto comp_gfs =
- std::make_shared<LGFS>(entity_set, finite_element_map);
+ auto comp_gfs = std::make_shared<LGFS>(entity_set, finite_element_map);
comp_gfs->name(name);
return comp_gfs;
@@ -149,8 +145,8 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
auto
-ModelDiffusionReaction<Traits>::
- setup_domain_grid_function_space(std::vector<std::string> comp_names) const
+ModelDiffusionReaction<Traits>::setup_domain_grid_function_space(
+ std::vector<std::string> comp_names) const
{
_logger.debug("setup domain grid function space"_fmt);
@@ -166,8 +162,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_grid_function_space()
+ModelDiffusionReaction<Traits>::setup_grid_function_space()
{
// get component names
auto& operator_splitting_config = _config.sub("operator");
@@ -201,8 +196,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_coefficient_vectors()
+ModelDiffusionReaction<Traits>::setup_coefficient_vectors()
{
_logger.debug("setup coefficient vector"_fmt);
for (auto& [op, state] : _states) {
@@ -217,8 +211,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_constraints()
+ModelDiffusionReaction<Traits>::setup_constraints()
{
_logger.debug("setup constraints"_fmt);
@@ -240,8 +233,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
auto
-ModelDiffusionReaction<Traits>::
- setup_local_operator(std::size_t i) const
+ModelDiffusionReaction<Traits>::setup_local_operator(std::size_t i) const
{
_logger.trace("setup local operators {}"_fmt, i);
@@ -260,8 +252,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_local_operators()
+ModelDiffusionReaction<Traits>::setup_local_operators()
{
_logger.trace("setup local operators"_fmt);
@@ -277,8 +268,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_grid_operators()
+ModelDiffusionReaction<Traits>::setup_grid_operators()
{
_logger.debug("setup grid operators"_fmt);
_spatial_grid_operators.clear();
@@ -343,8 +333,7 @@ ModelDiffusionReaction<Traits>::setup_solvers()
template<class Traits>
void
-ModelDiffusionReaction<Traits>::
- setup_vtk_writer()
+ModelDiffusionReaction<Traits>::setup_vtk_writer()
{
_logger.debug("setup vtk writer"_fmt);
@@ -370,8 +359,7 @@ ModelDiffusionReaction<Traits>::
template<class Traits>
void
-ModelDiffusionReaction<Traits>::setup(
- ModelSetupPolicy setup_policy)
+ModelDiffusionReaction<Traits>::setup(ModelSetupPolicy setup_policy)
{
if (setup_policy >= ModelSetupPolicy::GridFunctionSpace)
setup_grid_function_space();
@@ -391,8 +379,7 @@ ModelDiffusionReaction<Traits>::setup(
template<class Traits>
void
-ModelDiffusionReaction<Traits>::write_states()
- const
+ModelDiffusionReaction<Traits>::write_states() const
{
for (auto& [op, state] : _states) {
auto& x = state.coefficients;
@@ -408,8 +395,7 @@ ModelDiffusionReaction<Traits>::write_states()
using ET = decltype(etity_transformation);
using Predicate = PDELab::vtk::DefaultPredicate;
- using Data =
- PDELab::vtk::DGFTreeCommonData<GFS, X, Predicate, GV, ET>;
+ using Data = PDELab::vtk::DGFTreeCommonData<GFS, X, Predicate, GV, ET>;
std::shared_ptr<Data> data =
std::make_shared<Data>(*gfs, *x, _grid_view, etity_transformation);
PDELab::vtk::OutputCollector<SW, Data> collector(*_sequential_writer, data);
diff --git a/dune/copasi/model_diffusion_reaction.hh b/dune/copasi/model/diffusion_reaction.hh
similarity index 90%
rename from dune/copasi/model_diffusion_reaction.hh
rename to dune/copasi/model/diffusion_reaction.hh
index de8781d31ea2f4052a1533f1451e21436cab77fa..fd9ab700ec06f9d6e9d3ac8324db74b82a59d619 100644
--- a/dune/copasi/model_diffusion_reaction.hh
+++ b/dune/copasi/model/diffusion_reaction.hh
@@ -1,14 +1,13 @@
#ifndef DUNE_COPASI_MODEL_DIFFUSION_REACTION_HH
#define DUNE_COPASI_MODEL_DIFFUSION_REACTION_HH
-#include <dune/copasi/coefficient_mapper.hh>
+#include <dune/copasi/common/coefficient_mapper.hh>
+#include <dune/copasi/common/enum.hh>
#include <dune/copasi/concepts/grid.hh>
-#include <dune/copasi/enum.hh>
-#include <dune/copasi/grid_function_writer.hh>
-#include <dune/copasi/local_operator.hh>
-#include <dune/copasi/model_base.hh>
-#include <dune/copasi/model_state.hh>
-#include <dune/copasi/multidomain_local_finite_element_map.hh>
+#include <dune/copasi/finite_element/multidomain_local_finite_element_map.hh>
+#include <dune/copasi/model/base.hh>
+#include <dune/copasi/model/local_operator.hh>
+#include <dune/copasi/model/state.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/constraints/conforming.hh>
@@ -28,6 +27,15 @@
namespace Dune::Copasi {
+/**
+ * @brief Traits for diffusion reaction models
+ *
+ * @tparam G Grid
+ * @tparam GV Grid View
+ * @tparam FEMorder Order of the finite element method
+ * @tparam OT PDELab ordering tag
+ * @tparam JM Jacobian method
+ */
template<class G,
class GV = typename G::Traits::LeafGridView,
int FEMorder = 1,
@@ -44,9 +52,9 @@ struct ModelDiffusionReactionTraits
/**
* @brief Class for diffusion-reaction models.
+ * @todo Make this class work as stand-alone again
*
- * @tparam components Number of components
- * @tparam Param Parameterization class
+ * @tparam Traits Class that define static policies on the model
*/
template<class Traits>
class ModelDiffusionReaction : public ModelBase
@@ -209,16 +217,16 @@ public:
void step();
/**
- * @brief Get the model state
+ * @brief Get mutable model states
*
- * @return Model state
+ * @return Model states
*/
std::map<std::size_t, State> states() { return _states; }
/**
- * @brief Get the model state
+ * @brief Get constat model states
*
- * @return Model state
+ * @return Constant model states
*/
std::map<std::size_t, ConstState> const_states() const
{
@@ -228,9 +236,9 @@ public:
}
/**
- * @brief Get the model state
+ * @brief Get constat model states
*
- * @return Model state
+ * @return Constant model states
*/
std::map<std::size_t, ConstState> states() const { return const_states(); }
diff --git a/dune/copasi/local_operator.hh b/dune/copasi/model/local_operator.hh
similarity index 73%
rename from dune/copasi/local_operator.hh
rename to dune/copasi/model/local_operator.hh
index 3f43b1f118af1257a4211d00ea97510efbc62524..def904307167c96b4e002ece2d80935787b77aef 100644
--- a/dune/copasi/local_operator.hh
+++ b/dune/copasi/model/local_operator.hh
@@ -1,8 +1,9 @@
#ifndef DUNE_COPASI_LOCAL_OPERATOR_DIFFUSION_REACTION_HH
#define DUNE_COPASI_LOCAL_OPERATOR_DIFFUSION_REACTION_HH
-#include <dune/copasi/coefficient_mapper.hh>
-#include <dune/copasi/enum.hh>
+#include <dune/copasi/common/coefficient_mapper.hh>
+#include <dune/copasi/common/enum.hh>
+#include <dune/copasi/common/pdelab_expression_adapter.hh>
#include <dune/pdelab/common/quadraturerules.hh>
#include <dune/pdelab/localoperator/flags.hh>
@@ -16,12 +17,26 @@
#include <dune/common/fvector.hh>
-#include <dune/copasi/pdelab_expression_adapter.hh>
-
#include <set>
namespace Dune::Copasi {
+/**
+ * @brief This class describes a PDELab local operator for diffusion
+ * reaction systems.
+ * @details This class describre the operatrions for local integrals required
+ * for diffusion reaction system. The operator is only valid for
+ * entities contained in the grid view. The local finite element is
+ * used for caching shape function evaluations. The coefficient
+ * mapper is the interface to fetch date from either local or
+ * external coefficient vectors. And the jacobian method switches
+ * between numerical and analytical jacobians.
+ *
+ * @tparam GV Grid View
+ * @tparam LFE Local Finite Element
+ * @tparam CM Coefficient Mapper
+ * @tparam JM Jacobian Method
+ */
template<class GV,
class LFE,
class CM = DefaultCoefficientMapper,
@@ -106,14 +121,25 @@ public:
//! residual assembly flags
static constexpr bool doAlphaVolume = true;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @todo Make integration order variable depending on user requirements
+ * and polynomail order of the local finite element
+ *
+ * @param[in] grid_view The grid view where this local operator is
+ * valid
+ * @param[in] config The configuration tree
+ * @param[in] finite_element The local finite element
+ * @param[in] id_operator The index of this operator
+ */
LocalOperatorDiffusionReaction(GridView grid_view,
const ParameterTree& config,
const LocalFiniteElement& finite_element,
std::size_t id_operator)
: _basis_size(finite_element.localBasis().size())
, _components(config.sub("reaction").getValueKeys().size())
- , _rule(QuadratureRules<RF, dim>::rule(finite_element.type(),
- 3)) // TODO: make order variable
+ , _rule(QuadratureRules<RF, dim>::rule(finite_element.type(), 3))
, _diffusion_gf(_components)
, _reaction_gf(_components)
, _jacobian_gf(_components * _components)
@@ -215,12 +241,32 @@ public:
_logger.debug("LocalOperatorDiffusionReaction constructed"_fmt);
}
+ /**
+ * @brief Updates the coefficient mapper with given states.
+ *
+ * @param[in] states A map from operator index to states
+ *
+ * @tparam States Map from index to states
+ */
template<class States>
void update(const States& states)
{
_coefficient_mapper.update(states);
}
+ /**
+ * @brief Pattern volume
+ * @details This method links degrees of freedom between trial and test
+ * spaces taking into account the structure of the reaction term
+ *
+ * @param[in] lfsu The trial local function space
+ * @param[in] lfsv The test local function space
+ * @param pattern The local pattern
+ *
+ * @tparam LFSU The trial local function space
+ * @tparam LFSV The test local function space
+ * @tparam LocalPattern The local pattern
+ */
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume(const LFSU& lfsu,
const LFSV& lfsv,
@@ -238,6 +284,11 @@ public:
pattern.addLink(lfsv.child(i), k, lfsu.child(j), l);
}
+ /**
+ * @brief Sets the time.
+ *
+ * @param[in] t The new time
+ */
void setTime(double t)
{
Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>::setTime(t);
@@ -252,6 +303,25 @@ public:
gf->set_time(t);
}
+ /**
+ * @brief The jacobian volume integral for matrix free operations
+ * @details This only switches between the actual implementation (in
+ * _jacobian_apply_volume) and the numerical jacobian
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] z The local position in the trial space to which to apply
+ * the Jacobian.
+ * @param[in] lfsv The test local function space
+ * @param r The resulting vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The resulting vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
@@ -268,6 +338,57 @@ public:
}
}
+ /**
+ * @brief The jacobian volume integral for matrix free operations (linear
+ * variant)
+ * @details This only switches between the actual implementation (in
+ * _jacobian_apply_volume) and the numerical jacobian
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param r The resulting vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The resulting vector
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void jacobian_apply_volume(const EG& eg,
+ const LFSU& lfsu,
+ const X& x,
+ const LFSV& lfsv,
+ R& r) const
+ {
+ if constexpr (JM == JacobianMethod::Numerical) {
+ PDELab::NumericalJacobianApplyVolume<LocalOperatorDiffusionReaction>::
+ jacobian_apply_volume(eg, lfsu, x, lfsv, r);
+ return;
+ } else {
+ _jacobian_apply_volume(eg, lfsu, x, x, lfsv, r);
+ }
+ }
+
+ /**
+ * @brief The jacobian volume integral for matrix free operations
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] z The local position in the trial space to which to apply
+ * the Jacobian.
+ * @param[in] lfsv The test local function space
+ * @param r The resulting vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The resulting vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void _jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
@@ -361,7 +482,21 @@ public:
}
}
- //! jacobian contribution of volume term
+ /**
+ * @brief The jacobian volume integral
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param mat The local matrix
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam M The local matrix
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume(const EG& eg,
const LFSU& lfsu,
@@ -467,6 +602,21 @@ public:
}
}
+ /**
+ * @brief The volume integral
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param r The local residual vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The local residual vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(const EG& eg,
const LFSU& lfsu,
@@ -476,24 +626,21 @@ public:
{
_jacobian_apply_volume(eg, lfsu, x, x, lfsv, r);
}
-
- //! apply local jacobian of the volume term -> linear variant
- template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
- void jacobian_apply_volume(const EG& eg,
- const LFSU& lfsu,
- const X& x,
- const LFSV& lfsv,
- R& r) const
- {
- if constexpr (JM == JacobianMethod::Numerical) {
- PDELab::NumericalJacobianApplyVolume<LocalOperatorDiffusionReaction>::
- jacobian_apply_volume(eg, lfsu, x, lfsv, r);
- return;
- }
- _jacobian_apply_volume(eg, lfsu, x, x, lfsv, r);
- }
};
+/**
+ * @brief This class describes a PDELab local operator for temporal part
+ * for diffusion reaction systems.
+ * @details This class describre the operatrions for local integrals required
+ * for diffusion reaction system. The operator is only valid for
+ * entities contained in the grid view. The local finite element is
+ * used for caching shape function evaluations. And the jacobian
+ * method switches between numerical and analytical jacobians.
+ *
+ * @tparam GV Grid View
+ * @tparam LFE Local Finite Element
+ * @tparam JM Jacobian Method
+ */
template<class GV, class LFE, JacobianMethod JM = JacobianMethod::Analytical>
class TemporalLocalOperatorDiffusionReaction
: public Dune::PDELab::LocalOperatorDefaultFlags
@@ -560,6 +707,18 @@ public:
//! residual assembly flags
static constexpr bool doAlphaVolume = true;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @todo Make integration order variable depending on user requirements
+ * and polynomail order of the local finite element
+ *
+ * @param[in] grid_view The grid view where this local operator is
+ * valid
+ * @param[in] config The configuration tree
+ * @param[in] finite_element The local finite element
+ * @param[in] id_operator The index of this operator
+ */
TemporalLocalOperatorDiffusionReaction(
GridView grid_view,
const ParameterTree& config,
@@ -567,8 +726,7 @@ public:
std::size_t id_operator)
: _basis_size(finite_element.size())
, _components(config.sub("reaction").getValueKeys().size())
- , _rule(QuadratureRules<RF, dim>::rule(finite_element.type(),
- 3)) // TODO: make order variable
+ , _rule(QuadratureRules<RF, dim>::rule(finite_element.type(), 3))
, _logger(Logging::Logging::componentLogger(config, "model"))
{
auto& config_operator = config.sub("operator");
@@ -602,6 +760,21 @@ public:
_logger.debug("TemporalLocalOperatorDiffusionReaction constructed"_fmt);
}
+ /**
+ * @brief The volume integral
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param r The local residual vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The local residual vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(const EG& eg,
const LFSU& lfsu,
@@ -643,12 +816,27 @@ public:
}
}
- template<typename EG, typename LFSU, typename X, typename LFSV, typename Mat>
+ /**
+ * @brief The jacobian volume integral
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param mat The local matrix
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam M The local matrix
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume(const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
- Mat& mat) const
+ M& mat) const
{
if constexpr (JM == JacobianMethod::Numerical) {
PDELab::NumericalJacobianVolume<TemporalLocalOperatorDiffusionReaction>::
@@ -683,6 +871,25 @@ public:
}
}
+ /**
+ * @brief The jacobian volume integral for matrix free operations
+ * @details This only switches between the actual implementation (in
+ * alpha_volume) and the numerical jacobian
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] z The local position in the trial space to which to apply
+ * the Jacobian.
+ * @param[in] lfsv The test local function space
+ * @param r The resulting vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The resulting vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
@@ -704,6 +911,24 @@ public:
alpha_volume(eg, lfsu, z, lfsv, r);
}
+ /**
+ * @brief The jacobian volume integral for matrix free operations (linear
+ * variant)
+ * @details This only switches between the actual implementation (in
+ * alpha_volume) and the numerical jacobian
+ *
+ * @param[in] eg The entity
+ * @param[in] lfsu The trial local function space
+ * @param[in] x The local coefficient vector
+ * @param[in] lfsv The test local function space
+ * @param r The resulting vector
+ *
+ * @tparam EG The entity
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector type
+ * @tparam LFSV The test local function space
+ * @tparam R The resulting vector
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
diff --git a/dune/copasi/model_multidomain_diffusion_reaction.cc b/dune/copasi/model/multidomain_diffusion_reaction.cc
similarity index 91%
rename from dune/copasi/model_multidomain_diffusion_reaction.cc
rename to dune/copasi/model/multidomain_diffusion_reaction.cc
index 7b87ef283e81bab5c59835be0a29e38406c925a9..3ac990e4ecc53c338eff63e0e3303c68eeaa6073 100644
--- a/dune/copasi/model_multidomain_diffusion_reaction.cc
+++ b/dune/copasi/model/multidomain_diffusion_reaction.cc
@@ -7,18 +7,18 @@
* file but a header which has to be included when compiling.
*/
+#include <dune/copasi/common/muparser_data_handler.hh>
+#include <dune/copasi/model/multidomain_diffusion_reaction.hh>
+
#include <dune/pdelab/common/functionutilities.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionadapter.hh>
-#include <dune/copasi/model_multidomain_diffusion_reaction.hh>
-#include <dune/copasi/muparser_data_handler.hh>
-
namespace Dune::Copasi {
template<class Traits>
-ModelMultiDomainDiffusionReaction<Traits>::
- ModelMultiDomainDiffusionReaction(std::shared_ptr<Grid> grid,
- const Dune::ParameterTree& config)
+ModelMultiDomainDiffusionReaction<Traits>::ModelMultiDomainDiffusionReaction(
+ std::shared_ptr<Grid> grid,
+ const Dune::ParameterTree& config)
: ModelBase(config)
, _solver_logger(Logging::Logging::componentLogger(config, "solver"))
, _config(config)
@@ -92,16 +92,14 @@ ModelMultiDomainDiffusionReaction<Traits>::
}
template<class Traits>
-ModelMultiDomainDiffusionReaction<Traits>::
- ~ModelMultiDomainDiffusionReaction()
+ModelMultiDomainDiffusionReaction<Traits>::~ModelMultiDomainDiffusionReaction()
{
_logger.debug("ModelMultiDomainDiffusionReaction deconstructed"_fmt);
}
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_grid_function_spaces()
+ModelMultiDomainDiffusionReaction<Traits>::setup_grid_function_spaces()
{
_logger.debug("setup grid function space"_fmt);
@@ -155,8 +153,7 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_coefficient_vectors()
+ModelMultiDomainDiffusionReaction<Traits>::setup_coefficient_vectors()
{
for (auto& [op, state] : _states) {
_logger.debug("setup coefficient vector for operator {}"_fmt, op);
@@ -171,8 +168,7 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_constraints()
+ModelMultiDomainDiffusionReaction<Traits>::setup_constraints()
{
_logger.debug("setup constraints"_fmt);
@@ -198,8 +194,8 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
auto
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_local_operator(std::size_t i) const
+ModelMultiDomainDiffusionReaction<Traits>::setup_local_operator(
+ std::size_t i) const
{
_logger.trace("setup local operators {}"_fmt, i);
@@ -218,8 +214,7 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_local_operators()
+ModelMultiDomainDiffusionReaction<Traits>::setup_local_operators()
{
_logger.trace("setup local operators"_fmt);
@@ -235,8 +230,7 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_grid_operators()
+ModelMultiDomainDiffusionReaction<Traits>::setup_grid_operators()
{
_logger.debug("setup grid operators"_fmt);
_spatial_grid_operators.clear();
@@ -250,9 +244,9 @@ ModelMultiDomainDiffusionReaction<Traits>::
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());
+ max_comps = std::max(max_comps, gfs->child(i).degree());
- MBE mbe((int)pow(3, dim)*max_comps);
+ MBE mbe((int)pow(3, dim) * max_comps);
_logger.trace("create spatial grid operator {}"_fmt, op);
_spatial_grid_operators[op] = std::make_shared<GOS>(
@@ -299,8 +293,7 @@ ModelMultiDomainDiffusionReaction<Traits>::setup_solvers()
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- setup_vtk_writer()
+ModelMultiDomainDiffusionReaction<Traits>::setup_vtk_writer()
{
_logger.debug("setup vtk writer"_fmt);
const auto& compartments = _config.sub("compartments").getValueKeys();
@@ -343,16 +336,14 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- suggest_timestep(double dt)
+ModelMultiDomainDiffusionReaction<Traits>::suggest_timestep(double dt)
{
DUNE_THROW(NotImplemented, "not implemented");
}
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::setup(
- ModelSetupPolicy setup_policy)
+ModelMultiDomainDiffusionReaction<Traits>::setup(ModelSetupPolicy setup_policy)
{
_logger.trace("setup operator started"_fmt);
@@ -445,8 +436,8 @@ ModelMultiDomainDiffusionReaction<Traits>::step()
template<class Traits>
auto
-ModelMultiDomainDiffusionReaction<Traits>::
- get_data_handler(std::map<std::size_t, State> states) const
+ModelMultiDomainDiffusionReaction<Traits>::get_data_handler(
+ std::map<std::size_t, State> states) const
{
std::vector<std::map<std::size_t, std::shared_ptr<DataHandler>>> data(
_domains);
@@ -469,18 +460,17 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::
- update_data_handler()
+ModelMultiDomainDiffusionReaction<Traits>::update_data_handler()
{
_data = get_data_handler(_states);
}
template<class Traits>
auto
-ModelMultiDomainDiffusionReaction<Traits>::
- get_grid_function(const std::map<std::size_t, State>& states,
- std::size_t domain,
- std::size_t comp) const
+ModelMultiDomainDiffusionReaction<Traits>::get_grid_function(
+ const std::map<std::size_t, State>& states,
+ std::size_t domain,
+ std::size_t comp) const
{
auto data = get_data_handler(states);
const auto& compartments = _config.sub("compartments").getValueKeys();
@@ -503,8 +493,7 @@ ModelMultiDomainDiffusionReaction<Traits>::
template<class Traits>
void
-ModelMultiDomainDiffusionReaction<Traits>::write_states()
- const
+ModelMultiDomainDiffusionReaction<Traits>::write_states() const
{
const auto& compartments = _config.sub("compartments").getValueKeys();
diff --git a/dune/copasi/model_multidomain_diffusion_reaction.hh b/dune/copasi/model/multidomain_diffusion_reaction.hh
similarity index 77%
rename from dune/copasi/model_multidomain_diffusion_reaction.hh
rename to dune/copasi/model/multidomain_diffusion_reaction.hh
index 7f2f47465242e711a7598f3ade64e9fc4acf7e81..e95fc71b8a4526d4f8413a879b48b588319354be 100644
--- a/dune/copasi/model_multidomain_diffusion_reaction.hh
+++ b/dune/copasi/model/multidomain_diffusion_reaction.hh
@@ -1,13 +1,13 @@
#ifndef DUNE_COPASI_MODEL_MULTIDOMAIN_DIFFUSION_REACTION_HH
#define DUNE_COPASI_MODEL_MULTIDOMAIN_DIFFUSION_REACTION_HH
+#include <dune/copasi/common/enum.hh>
#include <dune/copasi/concepts/grid.hh>
-#include <dune/copasi/enum.hh>
-#include <dune/copasi/local_operator_multidomain.hh>
-#include <dune/copasi/model_base.hh>
-#include <dune/copasi/model_diffusion_reaction.cc>
-#include <dune/copasi/model_diffusion_reaction.hh>
-#include <dune/copasi/multidomain_entity_transformation.hh>
+#include <dune/copasi/grid/multidomain_entity_transformation.hh>
+#include <dune/copasi/model/base.hh>
+#include <dune/copasi/model/diffusion_reaction.cc>
+#include <dune/copasi/model/diffusion_reaction.hh>
+#include <dune/copasi/model/multidomain_local_operator.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/backend/istl/novlpistlsolverbackend.hh>
@@ -20,6 +20,14 @@
namespace Dune::Copasi {
+/**
+ * @brief Traits for diffusion reaction models in multigrid domains
+ *
+ * @tparam G Grid
+ * @tparam FEMorder Order of the finite element method
+ * @tparam OT PDELab ordering tag
+ * @tparam JM Jacobian method
+ */
template<class G,
int FEMorder = 1,
class OT = PDELab::EntityBlockedOrderingTag,
@@ -32,6 +40,11 @@ struct ModelMultiDomainDiffusionReactionTraits
static constexpr JacobianMethod jacobian_method = JM;
};
+/**
+ * @brief Class for diffusion-reaction models in multigrid domains.
+ *
+ * @tparam Traits Class that define static policies on the model
+ */
template<class Traits>
class ModelMultiDomainDiffusionReaction : public ModelBase
{
@@ -163,28 +176,44 @@ class ModelMultiDomainDiffusionReaction : public ModelBase
DGFTreeLeafFunction<ComponentLFS, DataHandler, SubDomainGridView>;
public:
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] grid The grid
+ * @param[in] config The configuration
+ */
ModelMultiDomainDiffusionReaction(std::shared_ptr<Grid> grid,
const Dune::ParameterTree& config);
+ /**
+ * @brief Destroys the object.
+ */
~ModelMultiDomainDiffusionReaction();
/**
- * @brief Get the model state
+ * @brief Get mutable model states
*
- * @return Model state
+ * @return Model states
*/
std::map<std::size_t, State> states() { return _states; }
/**
- * @brief Get the model state
+ * @brief Get constat model states
*
- * @return Model state
+ * @return Constant model states
*/
std::map<std::size_t, ConstState> const_states() const
{
return const_states(_states);
}
+ /**
+ * @brief Get constat model states
+ *
+ * @param[in] states A map to mutable states
+ *
+ * @return Constant model states
+ */
std::map<std::size_t, ConstState> const_states(
const std::map<std::size_t, State>& states) const
{
@@ -194,21 +223,44 @@ public:
}
/**
- * @brief Get the model state
+ * @brief Get the model states
*
- * @return Model state
+ * @return Model states
*/
std::map<std::size_t, ConstState> states() const { return const_states(); }
+ /**
+ * @brief Setup function
+ * @details This class sets up the model to have completely defined the
+ * requested feature in the policy
+ *
+ * @param[in] setup_policy The setup policy
+ */
void setup(ModelSetupPolicy setup_policy = ModelSetupPolicy::All);
+ /**
+ * @copydoc ModelBase::suggest_timestep
+ */
void suggest_timestep(double dt) override;
+ /**
+ * @copydoc ModelBase::step
+ */
void step() override;
- auto get_grid_function(const std::map<std::size_t, State>&,
- std::size_t,
- std::size_t) const;
+ /**
+ * @brief Gets a grid function for a given component and a sub domain.
+ * @todo Make this function operate on constant states
+ *
+ * @param[in] states The model states
+ * @param[in] domain The domain
+ * @param[in] comp The component
+ *
+ * @return The grid function.
+ */
+ auto get_grid_function(const std::map<std::size_t, State>& states,
+ std::size_t domain,
+ std::size_t comp) const;
protected:
void setup_grid_function_spaces();
diff --git a/dune/copasi/local_operator_multidomain.hh b/dune/copasi/model/multidomain_local_operator.hh
similarity index 80%
rename from dune/copasi/local_operator_multidomain.hh
rename to dune/copasi/model/multidomain_local_operator.hh
index 4f800788ee95d8da3d8dde9768bd4c1b02908f22..17208f68e3e61dc519deebe8bdad840fa9f373d4 100644
--- a/dune/copasi/local_operator_multidomain.hh
+++ b/dune/copasi/model/multidomain_local_operator.hh
@@ -1,9 +1,9 @@
#ifndef DUNE_COPASI_LOCAL_OPERATOR_MULTIDOMAIN_DIFFUSION_REACTION_HH
#define DUNE_COPASI_LOCAL_OPERATOR_MULTIDOMAIN_DIFFUSION_REACTION_HH
+#include <dune/copasi/common/enum.hh>
#include <dune/copasi/concepts/grid.hh>
-#include <dune/copasi/enum.hh>
-#include <dune/copasi/local_operator.hh>
+#include <dune/copasi/model/local_operator.hh>
#include <dune/pdelab/localoperator/numericaljacobian.hh>
#include <dune/pdelab/localoperator/numericaljacobianapply.hh>
@@ -15,6 +15,24 @@
namespace Dune::Copasi {
+/**
+ * @brief This class describes a PDELab local operator for multi domain
+ * diffusion reaction systems.
+ * @details This class describre the operatrions for local integrals required
+ * for diffusion reaction system. The operator is only valid for
+ * entities contained in the grid. The local finite element is used
+ * for caching shape function evaluations. The coefficient mapper is
+ * the interface to fetch date from either local or external
+ * coefficient vectors. And the jacobian method switches between
+ * numerical and analytical jacobians. This local operator creates
+ * internally an individual local operator for every subdomain in
+ * the grid
+ *
+ * @tparam Grid The grid
+ * @tparam LocalFiniteElement Local Finite Element
+ * @tparam CM Coefficient Mapper
+ * @tparam JM Jacobian Method
+ */
template<class Grid,
class LocalFiniteElement,
class CM = DefaultCoefficientMapper,
@@ -94,6 +112,14 @@ public:
//! residual assembly flags
static constexpr bool doAlphaSkeleton = true;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] grid The grid
+ * @param[in] config The configuration
+ * @param[in] finite_element The local finite element
+ * @param[in] id_operator The index of this operator
+ */
LocalOperatorMultiDomainDiffusionReaction(
std::shared_ptr<const Grid> grid,
const ParameterTree& config,
@@ -146,6 +172,13 @@ public:
}
}
+ /**
+ * @brief Updates the coefficient mapper with given states.
+ *
+ * @param[in] states A map from operator index to states
+ *
+ * @tparam States Map from index to states
+ */
template<class States>
void update(const States& states)
{
@@ -153,7 +186,9 @@ public:
_local_operator[i]->update(states);
}
- // define sparsity pattern of operator representation
+ /**
+ * @copydoc LocalOperatorDiffusionReaction::pattern_volume
+ */
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume(const LFSU& lfsu,
const LFSV& lfsv,
@@ -164,6 +199,23 @@ public:
}
}
+ /**
+ * @brief Pattern sckeleton
+ * @details This method links degrees of freedom between trial and test
+ * spaces at entities intersection taking into account the
+ * structure of the reaction term
+ *
+ * @param[in] lfsu_i The inside trial local function space
+ * @param[in] lfsv_i The inside test local function space
+ * @param[in] lfsu_o The outside trial local function space
+ * @param[in] lfsv_o The outside test local function space
+ * @param pattern_io The inside-outside local pattern
+ * @param pattern_oi The outside-inside local pattern
+ *
+ * @tparam LFSU The trial local function space
+ * @tparam LFSV The test local function space
+ * @tparam LocalPattern The local pattern
+ */
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_skeleton(const LFSU& lfsu_i,
const LFSV& lfsv_i,
@@ -235,6 +287,11 @@ public:
}
}
+ /**
+ * @brief Sets the time.
+ *
+ * @param[in] t The new time
+ */
void setTime(double t)
{
Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>::setTime(t);
@@ -242,6 +299,11 @@ public:
lp->setTime(t);
}
+ /**
+ * @copydoc LocalOperatorDiffusionReaction::jacobian_apply_volume
+ * @details This particular operator does a jacobian apply volume for the
+ * LocalOperatorDiffusionReaction corresponding to incoming entity
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
@@ -262,7 +324,34 @@ public:
}
}
- //! jacobian contribution of volume term
+ /**
+ * @copydoc LocalOperatorDiffusionReaction::jacobian_apply_volume
+ * @details This particular operator does a jacobian apply volume for the
+ * LocalOperatorDiffusionReaction corresponding to incoming entity
+ */
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void jacobian_apply_volume(const EG& eg,
+ const LFSU& lfsu,
+ const X& x,
+ const LFSV& lfsv,
+ R& r) const
+ {
+ for (std::size_t i = 0; i < _size; ++i) {
+ if (lfsu.child(i).size() > 0) {
+ const auto& sub_lfsu = lfsu.child(i);
+ const auto& sub_lfsv = lfsv.child(i);
+ if (sub_lfsu.size() == 0)
+ continue;
+ _local_operator[i]->jacobian_apply_volume(eg, sub_lfsu, x, sub_lfsv, r);
+ }
+ }
+ }
+
+ /**
+ * @copydoc LocalOperatorDiffusionReaction::jacobian_volume
+ * @details This particular operator does a jacobian volume for the
+ * LocalOperatorDiffusionReaction corresponding to incoming entity
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume(const EG& eg,
const LFSU& lfsu,
@@ -281,6 +370,11 @@ public:
}
}
+ /**
+ * @copydoc LocalOperatorDiffusionReaction::alpha_volume
+ * @details This particular operator does a alpha volume for the
+ * LocalOperatorDiffusionReaction corresponding to incoming entity
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(const EG& eg,
const LFSU& lfsu,
@@ -300,26 +394,28 @@ public:
}
}
- //! apply local jacobian of the volume term -> linear variant
- template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
- void jacobian_apply_volume(const EG& eg,
- const LFSU& lfsu,
- const X& x,
- const LFSV& lfsv,
- R& r) const
- {
- for (std::size_t i = 0; i < _size; ++i) {
- if (lfsu.child(i).size() > 0) {
- const auto& sub_lfsu = lfsu.child(i);
- const auto& sub_lfsv = lfsv.child(i);
- if (sub_lfsu.size() == 0)
- continue;
- _local_operator[i]->jacobian_apply_volume(eg, sub_lfsu, x, sub_lfsv, r);
- }
- }
- }
-
- // skeleton integral depending on test and ansatz functions
+ /**
+ * @brief The skeleton integral
+ * @details This integral is only performed at the interface between
+ * different domains. Currently it has the form of
+ * dichlet-dirichlet boundary condition between domains
+ *
+ * @param[in] ig The intersection
+ * @param[in] lfsu_i The inside trial local function space
+ * @param[in] x_i The inside local coefficient vector
+ * @param[in] lfsv_i The inside test local function space
+ * @param[in] lfsu_o The outside trial local function space
+ * @param[in] x_o The outside local coefficient vector
+ * @param[in] lfsv_o The outside test local function space
+ * @param r_i The inside residual vector
+ * @param r_o The outside residual vector
+ *
+ * @tparam IG The indersection
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector
+ * @tparam LFSV The test local function space
+ * @tparam R The residual vector
+ */
template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_skeleton(const IG& ig,
const LFSU& lfsu_i,
@@ -451,6 +547,28 @@ public:
}
}
+ /**
+ * @brief The jacobian skeleton integral
+ * @copydetails alpha_skeleton
+ *
+ * @param[in] ig The intersection
+ * @param[in] lfsu_i The inside trial local function space
+ * @param[in] x_i The inside local coefficient vector
+ * @param[in] lfsv_i The inside test local function space
+ * @param[in] lfsu_o The outside trial local function space
+ * @param[in] x_o The outside local coefficient vector
+ * @param[in] lfsv_o The outside test local function space
+ * @param mat_ii The local inside-inside matrix
+ * @param mat_io The local inside-outside matrix
+ * @param mat_oi The local outside-inside matrix
+ * @param mat_oo The local outside-outside matrix
+ *
+ * @tparam IG The indersection
+ * @tparam LFSU The trial local function space
+ * @tparam X The local coefficient vector
+ * @tparam LFSV The test local function space
+ * @tparam J The local jacobian matrix
+ */
template<typename IG, typename LFSU, typename X, typename LFSV, typename J>
void jacobian_skeleton(const IG& ig,
const LFSU& lfsu_i,
@@ -640,52 +758,23 @@ public:
}
}
}
-
- template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
- void jacobian_apply_skeleton(const IG& ig,
- const LFSU& lfsu_i,
- const X& x_i,
- const LFSV& lfsv_i,
- const LFSU& lfsu_o,
- const X& x_o,
- const LFSV& lfsv_o,
- Y& y_i,
- Y& y_o) const
- {
- if constexpr (JM == JacobianMethod::Numerical) {
- PDELab::NumericalJacobianApplySkeleton<
- LocalOperatorMultiDomainDiffusionReaction>::
- jacobian_apply_skeleton(
- ig, lfsu_i, x_i, lfsv_i, lfsu_o, x_o, lfsv_o, y_i, y_o);
- return;
- }
- DUNE_THROW(NotImplemented, "Analytic jacobian is not implemented");
- }
-
- template<typename IG, typename LFSU, typename X, typename LFSV, typename Y>
- void jacobian_apply_skeleton(const IG& ig,
- const LFSU& lfsu_i,
- const X& x_i,
- const X& z_i,
- const LFSV& lfsv_i,
- const LFSU& lfsu_o,
- const X& x_o,
- const X& z_o,
- const LFSV& lfsv_o,
- Y& y_i,
- Y& y_o) const
- {
- if constexpr (JM == JacobianMethod::Numerical) {
- PDELab::NumericalJacobianApplySkeleton<
- LocalOperatorMultiDomainDiffusionReaction>::
- jacobian_apply_skeleton(
- ig, lfsu_i, x_i, z_i, lfsv_i, lfsu_o, x_o, z_o, lfsv_o, y_i, y_o);
- return;
- }
- DUNE_THROW(NotImplemented, "Analytic jacobian is not implemented");
- }
};
+/**
+ * @brief This class describes a PDELab temporal local operator for multi
+ * domain diffusion reaction systems.
+ * @details This class describre the operatrions for local integrals required
+ * for diffusion reaction system. The operator is only valid for
+ * entities contained in the grid. The local finite element is used
+ * for caching shape function evaluations.And the jacobian method
+ * switches between numerical and analytical jacobians. This local
+ * operator creates internally an individual local operator for
+ * every subdomain in the grid
+ *
+ * @tparam Grid The grid
+ * @tparam LocalFiniteElement The local finite element
+ * @tparam JM The jacobian method
+ */
template<class Grid,
class LocalFiniteElement,
JacobianMethod JM = JacobianMethod::Analytical>
@@ -710,6 +799,14 @@ public:
//! residual assembly flags
static constexpr bool doAlphaVolume = true;
+ /**
+ * @brief Constructs a new instance.
+ *
+ * @param[in] grid The grid
+ * @param[in] config The configuration
+ * @param[in] finite_element The local finite element
+ * @param[in] id_operator The index of this operator
+ */
TemporalLocalOperatorMultiDomainDiffusionReaction(
std::shared_ptr<const Grid> grid,
const ParameterTree& config,
@@ -733,7 +830,9 @@ public:
}
}
- // define sparsity pattern of operator representation
+ /**
+ * @copydoc TemporalLocalOperatorDiffusionReaction::pattern_volume
+ */
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume(const LFSU& lfsu,
const LFSV& lfsv,
@@ -743,6 +842,9 @@ public:
_local_operator[i]->pattern_volume(lfsu.child(i), lfsv.child(i), pattern);
}
+ /**
+ * @copydoc LocalOperatorMultiDomainDiffusionReaction::alpha_volume
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(const EG& eg,
const LFSU& lfsu,
@@ -761,6 +863,9 @@ public:
}
}
+ /**
+ * @copydoc LocalOperatorMultiDomainDiffusionReaction::jacobian_volume
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename Mat>
void jacobian_volume(const EG& eg,
const LFSU& lfsu,
@@ -779,6 +884,9 @@ public:
}
}
+ /**
+ * @copydoc LocalOperatorMultiDomainDiffusionReaction::jacobian_apply_volume
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
@@ -799,6 +907,9 @@ public:
}
}
+ /**
+ * @copydoc LocalOperatorMultiDomainDiffusionReaction::jacobian_apply_volume
+ */
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void jacobian_apply_volume(const EG& eg,
const LFSU& lfsu,
diff --git a/dune/copasi/model_state.hh b/dune/copasi/model/state.hh
similarity index 91%
rename from dune/copasi/model_state.hh
rename to dune/copasi/model/state.hh
index 94c3da8b71bcb9c46089b3c5dac6aa4968dfc1ee..131fca52099f4bb594c85e7ea471a3845d506fa7 100644
--- a/dune/copasi/model_state.hh
+++ b/dune/copasi/model/state.hh
@@ -7,9 +7,17 @@
namespace Dune::Copasi {
+/// forward declaration of model state
template<class G, class GFS, class X>
struct ModelState;
+/**
+ * @brief Constant model state container
+ *
+ * @tparam G Grid type
+ * @tparam GFS Grid function space type
+ * @tparam X Coefficients vector
+ */
template<class G, class GFS, class X>
struct ConstModelState
{
diff --git a/dune/copasi/muparser_data_handler.hh b/dune/copasi/muparser_data_handler.hh
deleted file mode 100644
index 7f60f1d877ad79efd50022ab4554dc3c5e0cd141..0000000000000000000000000000000000000000
--- a/dune/copasi/muparser_data_handler.hh
+++ /dev/null
@@ -1,67 +0,0 @@
-#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
diff --git a/dune/copasi/pdelab_callable_adapter.hh b/dune/copasi/pdelab_callable_adapter.hh
deleted file mode 100644
index e095af244207679336a86aa27f0fc028926a4e98..0000000000000000000000000000000000000000
--- a/dune/copasi/pdelab_callable_adapter.hh
+++ /dev/null
@@ -1,142 +0,0 @@
-#ifndef DUNE_COPASI_GRID_FUNCTION_CALLABLE_ADAPTER_HH
-#define DUNE_COPASI_GRID_FUNCTION_CALLABLE_ADAPTER_HH
-
-#include <dune/copasi/concepts/pdelab.hh>
-
-#include <dune/pdelab/common/function.hh>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/typetraits.hh>
-
-#include <type_traits>
-#include <utility>
-
-namespace Dune::Copasi {
-
-namespace {
-// namespaces unique to this translation unit/file
-
-template<class GV, class R>
-class GridFunctionTraitsFromRange;
-
-template<class GV, class RF, int dim>
-class GridFunctionTraitsFromRange<GV, Dune::FieldVector<RF, dim>>
- : public PDELab::GridFunctionTraits<GV, RF, dim, Dune::FieldVector<RF, dim>>
-{};
-
-template<class GV, class RF>
-class GridFunctionTraitsFromRange<GV, Dune::DynamicVector<RF>>
- : public PDELab::GridFunctionTraits<GV, RF, -1, Dune::DynamicVector<RF>>
-{};
-
-} // namespace
-
-template<typename GV, typename R, typename F>
-class GlobalCallableToGridFunctionAdapter
- : public Dune::PDELab::GridFunctionBase<
- GridFunctionTraitsFromRange<GV, R>,
- GlobalCallableToGridFunctionAdapter<GV, R, F>>
-{
-public:
- typedef GridFunctionTraitsFromRange<GV, R> Traits;
-
- //! construct from grid view
- GlobalCallableToGridFunctionAdapter(const GV& grid_view, const F& callable)
- : _gv(grid_view)
- , _f(callable)
- {}
-
- //! get a reference to the grid view
- inline const GV& getGridView() const { return _gv; }
-
- //! evaluate extended function on element
- inline void evaluate(const typename Traits::ElementType& e,
- const typename Traits::DomainType& xl,
- typename Traits::RangeType& y) const
- {
- typename Traits::DomainType xg = e.geometry().global(xl);
- y = _f(xg);
- }
-
-private:
- GV _gv;
- F _f;
-};
-
-template<typename GV, typename R, typename F>
-class LocalCallableToGridFunctionAdapter
- : public Dune::PDELab::GridFunctionBase<
- GridFunctionTraitsFromRange<GV, R>,
- LocalCallableToGridFunctionAdapter<GV, R, F>>
-{
-public:
- typedef GridFunctionTraitsFromRange<GV, R> Traits;
-
- //! construct from grid view
- LocalCallableToGridFunctionAdapter(const GV& grid_view, const F& callable)
- : _gv(grid_view)
- , _f(callable)
- {}
-
- //! get a reference to the grid view
- inline const GV& getGridView() const { return _gv; }
-
- //! evaluate extended function on element
- inline void evaluate(const typename Traits::ElementType& e,
- const typename Traits::DomainType& xl,
- typename Traits::RangeType& y) const
- {
- y = _f(e, xl);
- }
-
-private:
- GV _gv;
- F _f;
-};
-
-/** \brief Create PDELab GridFunction from a callable f(x) that expects a global
- * coordinate x */
-template<typename GV, typename F>
-auto
-makeGridFunctionFromCallable(const GV& gv, const F& f) ->
- typename std::enable_if<
- Concept::isPDELabGlobalCallable<GV, F>(),
- GlobalCallableToGridFunctionAdapter<
- GV,
- decltype(f(std::declval<typename GV::template Codim<
- 0>::Entity::Geometry::GlobalCoordinate>())),
- F>>::type
-{
- typedef typename GV::template Codim<0>::Entity::Geometry::GlobalCoordinate X;
- X x;
- typedef decltype(f(x)) Range;
- typedef GlobalCallableToGridFunctionAdapter<GV, Range, F> TheType;
- return TheType(gv, f);
-}
-
-/** \brief Create PDELab GridFunction from a callable f(e,x) that expects
- an entity e and a local coordinate x */
-template<typename GV, typename F>
-auto
-makeGridFunctionFromCallable(const GV& gv, const F& f) ->
- typename std::enable_if<
- Concept::isPDELabLocalCallable<GV, F>(),
- LocalCallableToGridFunctionAdapter<
- GV,
- decltype(f(std::declval<typename GV::template Codim<0>::Entity>(),
- std::declval<typename GV::template Codim<
- 0>::Entity::Geometry::LocalCoordinate>())),
- F>>::type
-{
- typedef typename GV::template Codim<0>::Entity E;
- E e;
- typedef typename E::Geometry::LocalCoordinate X;
- X x;
- typedef decltype(f(e, x)) Range;
- typedef LocalCallableToGridFunctionAdapter<GV, Range, F> TheType;
- return TheType(gv, f);
-}
-
-} // namespace Dune::Copasi
-
-#endif // DUNE_COPASI_GRID_FUNCTION_CALLABLE_ADAPTER_HH
\ No newline at end of file
diff --git a/dune/copasi/tiff_grayscale.hh b/dune/copasi/tiff_grayscale.hh
deleted file mode 100644
index 63c2639dea5ecc2a052818848a5d74fdfb499917..0000000000000000000000000000000000000000
--- a/dune/copasi/tiff_grayscale.hh
+++ /dev/null
@@ -1,158 +0,0 @@
-#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
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 9825ccfd9e6ce151b73f72069f174a411d4d60ff..e5d249db320f93243105d2fd4a3b94abc20b8754 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -1,5 +1,5 @@
-add_library(dune_copasi_lib STATIC ../dune/copasi/model_base.cc
- model_multidomain_diffusion_reaction.cc)
+add_library(dune_copasi_lib STATIC ../dune/copasi/model/base.cc
+ multidomain_diffusion_reaction.cc)
target_link_libraries(dune_copasi_lib PUBLIC ${DUNE_LIBS} TIFF::TIFF muparser::muparser)
diff --git a/lib/model_multidomain_diffusion_reaction.cc b/lib/model_multidomain_diffusion_reaction.cc
deleted file mode 100644
index b3fb2ebe95460a926c50fae17b3d6c1f33e818d2..0000000000000000000000000000000000000000
--- a/lib/model_multidomain_diffusion_reaction.cc
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/copasi/model_multidomain_diffusion_reaction.hh>
-#include <dune/copasi/model_multidomain_diffusion_reaction.cc>
-
-#include <dune/grid/multidomaingrid.hh>
-
-#include <dune/grid/uggrid.hh>
-
-namespace Dune{
-namespace Copasi{
-
-constexpr int dim = 2;
-using HostGrid = Dune::UGGrid<dim>;
-using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim,1>;
-using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid,MDGTraits>;
-
-using ModelTraits1 = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,1>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits1>;
-
-using ModelTraits2 = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,2>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits2>;
-
-using Ordering = Dune::PDELab::EntityBlockedOrderingTag;
-constexpr Dune::Copasi::JacobianMethod Jac = Dune::Copasi::JacobianMethod::Numerical;
-using ModelTraits1Num = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,1,Ordering,Jac>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits1Num>;
-
-} // namespace Dorie
-} // namespace Dune
diff --git a/lib/multidomain_diffusion_reaction.cc b/lib/multidomain_diffusion_reaction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..050efaf122f5eab0736059647af84a8a594a22c6
--- /dev/null
+++ b/lib/multidomain_diffusion_reaction.cc
@@ -0,0 +1,36 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/copasi/model/multidomain_diffusion_reaction.cc>
+#include <dune/copasi/model/multidomain_diffusion_reaction.hh>
+
+#include <dune/grid/multidomaingrid.hh>
+
+#include <dune/grid/uggrid.hh>
+
+namespace Dune {
+namespace Copasi {
+
+constexpr int dim = 2;
+using HostGrid = Dune::UGGrid<dim>;
+using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
+using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+
+using ModelTraits1 =
+ Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid, 1>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits1>;
+
+using ModelTraits2 =
+ Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid, 2>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2>;
+
+using Ordering = Dune::PDELab::EntityBlockedOrderingTag;
+constexpr Dune::Copasi::JacobianMethod Jac =
+ Dune::Copasi::JacobianMethod::Numerical;
+using ModelTraits1Num =
+ Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid, 1, Ordering, Jac>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits1Num>;
+
+} // namespace Dorie
+} // namespace Dune
diff --git a/src/dune_copasi.cc b/src/dune_copasi.cc
index e8068f4d83188ec2750116d3d25573ad0ca4e2d4..e52e7327648e3e441bf5d7658db34591735f12b9 100644
--- a/src/dune_copasi.cc
+++ b/src/dune_copasi.cc
@@ -1,35 +1,36 @@
#ifdef HAVE_CONFIG_H
-# include "config.h"
+#include "config.h"
#endif
-#include <dune/copasi/gmsh_reader.hh>
-#include <dune/copasi/model_diffusion_reaction.hh>
-#include <dune/copasi/model_diffusion_reaction.cc>
-#include <dune/copasi/model_multidomain_diffusion_reaction.hh>
-#include <dune/copasi/enum.hh>
+#include <dune/copasi/common/enum.hh>
+#include <dune/copasi/grid/multidomain_gmsh_reader.hh>
+#include <dune/copasi/model/diffusion_reaction.cc>
+#include <dune/copasi/model/diffusion_reaction.hh>
+#include <dune/copasi/model/multidomain_diffusion_reaction.hh>
-#include <dune/grid/multidomaingrid.hh>
#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/multidomaingrid.hh>
#include <dune/logging/logging.hh>
#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
-#include <dune/common/parallel/mpihelper.hh>
#include <iostream>
-int main(int argc, char** argv)
+int
+main(int argc, char** argv)
{
- try{
+ try {
// initialize mpi
auto& mpi_helper = Dune::MPIHelper::instance(argc, argv);
auto comm = mpi_helper.getCollectiveCommunication();
// Read and parse ini file
- if (argc!=2)
+ if (argc != 2)
DUNE_THROW(Dune::IOError, "Wrong number of arguments");
const std::string config_filename = argv[1];
@@ -38,7 +39,7 @@ int main(int argc, char** argv)
ptreeparser.readINITree(config_filename, config);
// initialize loggers
- Dune::Logging::Logging::init(comm,config.sub("logging"));
+ Dune::Logging::Logging::init(comm, config.sub("logging"));
using namespace Dune::Literals;
auto log = Dune::Logging::Logging::logger(config);
@@ -47,49 +48,49 @@ int main(int argc, char** argv)
// create a grid
constexpr int dim = 2;
using HostGrid = Dune::UGGrid<dim>;
- using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim,1>;
- using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid,MDGTraits>;
+ using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
+ using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
auto& grid_config = config.sub("grid");
- auto level = grid_config.get<int>("initial_level",0);
+ auto level = grid_config.get<int>("initial_level", 0);
log.info("Creating a rectangular grid in {}D"_fmt, dim);
auto grid_file = grid_config.get<std::string>("file");
- auto [grid_ptr,host_grid_ptr] = Dune::Copasi::GmshReader<Grid>::read(grid_file,config);
+ auto [grid_ptr, host_grid_ptr] =
+ Dune::Copasi::MultiDomainGmshReader<Grid>::read(grid_file, config);
log.debug("Applying global refinement of level: {}"_fmt, level);
grid_ptr->globalRefine(level);
auto& model_config = config.sub("model");
int order = model_config.get<int>("order");
-
- if (order == 1)
- {
+
+ if (order == 1) {
constexpr int Order = 1;
- using ModelTraits = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,Order>;
- Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(grid_ptr,model_config);
+ using ModelTraits =
+ Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid, Order>;
+ Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
+ grid_ptr, model_config);
model.run();
- }
- else if (order == 2)
- {
+ } else if (order == 2) {
constexpr int Order = 2;
- using ModelTraits = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,Order>;
- Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(grid_ptr,model_config);
+ using ModelTraits =
+ Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid, Order>;
+ Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
+ grid_ptr, model_config);
model.run();
- }
- else
- {
- DUNE_THROW(Dune::IOError,"Finite element order " << order << " is not supported by dune-copasi");
+ } else {
+ DUNE_THROW(Dune::IOError,
+ "Finite element order " << order
+ << " is not supported by dune-copasi");
}
return 0;
- }
- catch (Dune::Exception &e){
+ } catch (Dune::Exception& e) {
std::cerr << "Dune reported error: " << e << std::endl;
- }
- catch (...){
+ } catch (...) {
std::cerr << "Unknown exception thrown!" << std::endl;
}
}
diff --git a/test/test_concepts_grid.cc b/test/test_concepts_grid.cc
index f2e193ba8a3881c7b4009cc1a89348a02961f320..ad9f46d66b78282202d697aa1d076a1b3882f7e8 100644
--- a/test/test_concepts_grid.cc
+++ b/test/test_concepts_grid.cc
@@ -28,11 +28,14 @@ test_multidomaingrid()
// check that it is a multidomain grid with its respective subdomain grids
passed &= isGrid<Grid>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isMultiDomainGrid<Grid>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isSubDomainGrid<typename Grid::SubDomainGrid>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
if constexpr (i > 0)
// test for nested multidomains
@@ -44,7 +47,7 @@ test_multidomaingrid()
int
main(int argc, char** argv)
{
- [[maybe_unused]] auto& mpi_helper = Dune::MPIHelper::instance(argc, argv);
+ [[maybe_unused]] auto& mpi_helper = Dune::MPIHelper::instance(argc, argv);
try {
bool passed = true;
@@ -53,53 +56,73 @@ main(int argc, char** argv)
// check grids yaspgrid dim 1
passed &= isGrid<Dune::YaspGrid<1>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isMultiDomainGrid<Dune::YaspGrid<1>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isSubDomainGrid<Dune::YaspGrid<1>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= test_multidomaingrid<Dune::YaspGrid<1>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check grids yaspgrid dim 2
passed &= isGrid<Dune::YaspGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isMultiDomainGrid<Dune::YaspGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isSubDomainGrid<Dune::YaspGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= test_multidomaingrid<Dune::YaspGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check grids yaspgrid dim 5
passed &= isGrid<Dune::YaspGrid<5>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isMultiDomainGrid<Dune::YaspGrid<5>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isSubDomainGrid<Dune::YaspGrid<5>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= test_multidomaingrid<Dune::YaspGrid<5>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check grids uggrid dim 2
passed &= isGrid<Dune::UGGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isMultiDomainGrid<Dune::UGGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isSubDomainGrid<Dune::UGGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= test_multidomaingrid<Dune::UGGrid<2>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check grids uggrid dim 3
passed &= isGrid<Dune::UGGrid<3>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isMultiDomainGrid<Dune::UGGrid<3>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isSubDomainGrid<Dune::UGGrid<3>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= test_multidomaingrid<Dune::UGGrid<3>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
return not passed;
} catch (Dune::Exception& e) {
diff --git a/test/test_concepts_pdelab.cc b/test/test_concepts_pdelab.cc
index 394c7dbdcc7b93a62ca4042f9a99360ca9871e66..92ca4de05c32caec5619fb9e31442959fcbf5f27 100644
--- a/test/test_concepts_pdelab.cc
+++ b/test/test_concepts_pdelab.cc
@@ -9,8 +9,8 @@
#include <dune/grid/yaspgrid.hh>
-#include <dune/common/math.hh>
#include <dune/common/exceptions.hh>
+#include <dune/common/math.hh>
#include <cassert>
#include <complex>
@@ -47,18 +47,24 @@ main(int argc, char** argv)
// check functions
passed &= isPDELabFunction<F<double>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isPDELabFunction<F<int>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isPDELabFunction<F<std::complex<double>>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabFunction<double>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabFunction<int>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabFunction<std::complex<double>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// we need a grid for the following checks
int constexpr dimDomain = 3;
@@ -78,46 +84,58 @@ main(int argc, char** argv)
using PDELabLocalCallable = decltype(pdelab_local_callable);
passed &= isPDELabLocalCallable<GridView, PDELabLocalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isPDELabCallable<GridView, PDELabLocalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabLocalCallable<GridView, F<double>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabLocalCallable<GridView, double>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check global callable
auto pdelab_global_callable = [](const auto& x) { return x; };
using PDELabGlobalCallable = decltype(pdelab_global_callable);
passed &= isPDELabGlobalCallable<GridView, PDELabGlobalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= isPDELabCallable<GridView, PDELabGlobalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabGlobalCallable<GridView, F<double>>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabGlobalCallable<GridView, double>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabGlobalCallable<GridView, PDELabLocalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
passed &= not isPDELabLocalCallable<GridView, PDELabGlobalCallable>();
- if (not passed) DUNE_THROW(Dune::Exception, "");
+ if (not passed)
+ DUNE_THROW(Dune::Exception, "");
// check grid function from the pdelab creator from callables
- [[maybe_unused]] auto gf_from_local = Dune::PDELab::makeGridFunctionFromCallable(
- grid_view, pdelab_local_callable);
- [[maybe_unused]] auto gf_from_global = Dune::PDELab::makeGridFunctionFromCallable(
- grid_view, pdelab_global_callable);
+ [[maybe_unused]] auto gf_from_local =
+ Dune::PDELab::makeGridFunctionFromCallable(grid_view,
+ pdelab_local_callable);
+ [[maybe_unused]] auto gf_from_global =
+ Dune::PDELab::makeGridFunctionFromCallable(grid_view,
+ pdelab_global_callable);
// TODO: For some reason, they are not recoginized as grid functions
// passed &= isPDELabGridFunction<decltype(gf_from_local)>();
// if (not passed) DUNE_THROW(Dune::Exception, "");
// passed &= isPDELabGridFunction<decltype(gf_from_global)>();
// if (not passed) DUNE_THROW(Dune::Exception, "");
-
+
return not passed;
} catch (Dune::Exception& e) {
std::cerr << "Dune reported error: " << e << std::endl;
diff --git a/test/test_dynamic_local_finite_element.cc b/test/test_dynamic_local_finite_element.cc
index a78e47dcf1811eed4b18862d87d74edd8fa0a375..8553a7a8b48512d8dfcb645632b8581a372f72de 100644
--- a/test/test_dynamic_local_finite_element.cc
+++ b/test/test_dynamic_local_finite_element.cc
@@ -2,9 +2,9 @@
#include "config.h"
#endif
-#include <dune/copasi/dynamic_local_basis.hh>
-#include <dune/copasi/dynamic_local_coefficients.hh>
-#include <dune/copasi/dynamic_local_finite_element.hh>
+#include <dune/copasi/finite_element/dynamic_local_basis.hh>
+#include <dune/copasi/finite_element/dynamic_local_coefficients.hh>
+#include <dune/copasi/finite_element/dynamic_local_finite_element.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/lagrange.hh>
diff --git a/test/test_jacobian.cc b/test/test_jacobian.cc
index be667b3c6ebf518cf61e2944b86d6c7c4d1a282a..f40c1ba99db754e03466b1aebb9abc97222c2447 100644
--- a/test/test_jacobian.cc
+++ b/test/test_jacobian.cc
@@ -1,35 +1,36 @@
#ifdef HAVE_CONFIG_H
-# include "config.h"
+#include "config.h"
#endif
-#include <dune/copasi/gmsh_reader.hh>
-#include <dune/copasi/model_diffusion_reaction.hh>
-#include <dune/copasi/model_diffusion_reaction.cc>
-#include <dune/copasi/model_multidomain_diffusion_reaction.hh>
-#include <dune/copasi/enum.hh>
+#include <dune/copasi/common/enum.hh>
+#include <dune/copasi/grid/multidomain_gmsh_reader.hh>
+#include <dune/copasi/model/diffusion_reaction.cc>
+#include <dune/copasi/model/diffusion_reaction.hh>
+#include <dune/copasi/model/multidomain_diffusion_reaction.hh>
-#include <dune/grid/multidomaingrid.hh>
#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/multidomaingrid.hh>
#include <dune/logging/logging.hh>
#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
-#include <dune/common/parallel/mpihelper.hh>
#include <iostream>
-int main(int argc, char** argv)
+int
+main(int argc, char** argv)
{
- try{
+ try {
// initialize mpi
auto& mpi_helper = Dune::MPIHelper::instance(argc, argv);
auto comm = mpi_helper.getCollectiveCommunication();
// Read and parse ini file
- if (argc!=2)
+ if (argc != 2)
DUNE_THROW(Dune::IOError, "Wrong number of arguments");
const std::string config_filename = argv[1];
@@ -38,7 +39,7 @@ int main(int argc, char** argv)
ptreeparser.readINITree(config_filename, config);
// initialize loggers
- Dune::Logging::Logging::init(comm,config.sub("logging"));
+ Dune::Logging::Logging::init(comm, config.sub("logging"));
using namespace Dune::Literals;
auto log = Dune::Logging::Logging::logger(config);
@@ -47,17 +48,18 @@ int main(int argc, char** argv)
// create a grid
constexpr int dim = 2;
using HostGrid = Dune::UGGrid<dim>;
- using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim,1>;
- using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid,MDGTraits>;
+ using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
+ using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
auto& grid_config = config.sub("grid");
- auto level = grid_config.get<int>("initial_level",0);
+ auto level = grid_config.get<int>("initial_level", 0);
log.info("Creating a rectangular grid in {}D"_fmt, dim);
auto grid_file = grid_config.get<std::string>("file");
- auto [grid_ptr,host_grid_ptr] = Dune::Copasi::GmshReader<Grid>::read(grid_file,config);
+ auto [grid_ptr, host_grid_ptr] =
+ Dune::Copasi::MultiDomainGmshReader<Grid>::read(grid_file, config);
log.debug("Applying global refinement of level: {}"_fmt, level);
grid_ptr->globalRefine(level);
@@ -65,41 +67,41 @@ int main(int argc, char** argv)
auto& model_config = config.sub("model");
int order = model_config.get<int>("order");
std::string jacobian = model_config.get<std::string>("jacobian_type");
-
- if (order == 1)
- {
+
+ if (order == 1) {
using Ordering = Dune::PDELab::EntityBlockedOrderingTag;
constexpr int Order = 1;
- if (jacobian == "analytical")
- {
- constexpr Dune::Copasi::JacobianMethod Jac = Dune::Copasi::JacobianMethod::Analytical;
- using ModelTraits = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,Order,Ordering,Jac>;
- Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(grid_ptr,model_config);
+ if (jacobian == "analytical") {
+ constexpr Dune::Copasi::JacobianMethod Jac =
+ Dune::Copasi::JacobianMethod::Analytical;
+ using ModelTraits = Dune::Copasi::
+ ModelMultiDomainDiffusionReactionTraits<Grid, Order, Ordering, Jac>;
+ Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
+ grid_ptr, model_config);
model.run();
- }
- else if (jacobian == "numerical")
- {
- constexpr Dune::Copasi::JacobianMethod Jac = Dune::Copasi::JacobianMethod::Numerical;
- using ModelTraits = Dune::Copasi::ModelMultiDomainDiffusionReactionTraits<Grid,Order,Ordering,Jac>;
- Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(grid_ptr,model_config);
+ } else if (jacobian == "numerical") {
+ constexpr Dune::Copasi::JacobianMethod Jac =
+ Dune::Copasi::JacobianMethod::Numerical;
+ using ModelTraits = Dune::Copasi::
+ ModelMultiDomainDiffusionReactionTraits<Grid, Order, Ordering, Jac>;
+ Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
+ grid_ptr, model_config);
model.run();
+ } else {
+ DUNE_THROW(Dune::IOError,
+ "Jacobian type " << jacobian
+ << " is not supported by dune-copasi");
}
- else
- {
- DUNE_THROW(Dune::IOError,"Jacobian type " << jacobian << " is not supported by dune-copasi");
- }
- }
- else
- {
- DUNE_THROW(Dune::IOError,"Finite element order " << order << " is not supported by dune-copasi");
+ } else {
+ DUNE_THROW(Dune::IOError,
+ "Finite element order " << order
+ << " is not supported by dune-copasi");
}
return 0;
- }
- catch (Dune::Exception &e){
+ } catch (Dune::Exception& e) {
std::cerr << "Dune reported error: " << e << std::endl;
- }
- catch (...){
+ } catch (...) {
std::cerr << "Unknown exception thrown!" << std::endl;
}
}
diff --git a/test/test_parameter_tree_check.cc b/test/test_parameter_tree_check.cc
new file mode 100644
index 0000000000000000000000000000000000000000..bd4c7f8d79b3df2bf2513e672d848e6616c594fd
--- /dev/null
+++ b/test/test_parameter_tree_check.cc
@@ -0,0 +1,72 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+
+#include <dune/copasi/parameter_tree_check.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parametertreeparser.hh>
+
+int
+main(int argc, char** argv)
+{
+ bool failed = false;
+
+ try {
+
+ Dune::ParameterTree param_a, param_b, param_c;
+
+ param_a["a"] = "1";
+ param_a["A.a"] = "11";
+ param_b["b"] = "2";
+ param_c["C.c.a"] = "333";
+
+ if (not Dune::Copasi::eq(param_a,param_a))
+ DUNE_THROW(Dune::RangeError,"eq operator is wrong");
+
+ if (Dune::Copasi::eq(param_a,param_b))
+ DUNE_THROW(Dune::RangeError,"eq operator is wrong");
+
+ Dune::Copasi::add(param_a,param_b);
+
+ param_b["a"] = "1";
+ if (Dune::Copasi::eq(param_a,param_b))
+ DUNE_THROW(Dune::RangeError,"eq operator is wrong");
+
+ param_b["A.a"] = "11";
+ if (not Dune::Copasi::eq(param_a,param_b))
+ DUNE_THROW(Dune::RangeError,"eq operator is wrong");
+
+ Dune::Copasi::diff(param_a,param_b);
+ if (not Dune::Copasi::eq(param_a,{}))
+ DUNE_THROW(Dune::RangeError,"diff operator is wrong");
+
+ Dune::ParameterTree param_c_copy(param_c);
+ Dune::Copasi::add(param_c,param_b);
+ Dune::Copasi::add(param_a,param_b);
+ Dune::Copasi::diff(param_c,param_a);
+
+ if (not Dune::Copasi::eq(param_c,param_c_copy))
+ DUNE_THROW(Dune::RangeError,"diff operator is wrong");
+
+ const std::string config_filename = argv[1];
+ Dune::ParameterTree config;
+ Dune::ParameterTreeParser ptreeparser;
+ ptreeparser.readINITree(config_filename, config);
+
+ auto config_model = Dune::Copasi::get_model(config);
+ std::cout << "*******************Valid configuration model*******************" << std::endl;
+ config_model.report();
+
+ std::cout << "***************************Diff***************************" << std::endl;
+ Dune::Copasi::diff(config,config_model);
+ config.report();
+
+ return failed;
+ } catch (Dune::Exception& e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ } catch (...) {
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}
diff --git a/test/test_tiff_grayscale.cc b/test/test_tiff_grayscale.cc
index 2dbd9fe5a30703be4bdca93a1a37a2035a8bab8e..c3e85243173f70f645f488ac34304807ba2e766c 100644
--- a/test/test_tiff_grayscale.cc
+++ b/test/test_tiff_grayscale.cc
@@ -2,7 +2,7 @@
#include "config.h"
#endif
-#include <dune/copasi/tiff_grayscale.hh>
+#include <dune/copasi/common/tiff_grayscale.hh>
#include <dune/grid/yaspgrid.hh>