1 Introduction
Various simulation problems are posed on domains that are not open subsets of a Euclidean space. Frequently, such domains are surfaces or curves embedded in a higherdimensional Euclidean space. Equations on such domains, sometimes called geometric partial differential equations, comprise diffusion and transport on the surface
(Dziuk and Elliott, 2007b), flow problems (Nitschke et al., 2012; Reuther and Voigt, 2015), and phasefield equations (T. Witkowski, 2012). Sometimes, movement of the surface itself is modeled (Dziuk and Elliott, 2007a), and this movement may couple with processes on the surface (Gross and Reusken, 2011).As an additional difficulty, some boundary value problems are posed on domains that do not even have the structure of a topological manifold. That is, not every point of has a neighborhood that is homeomorphic to an open subset of . Figure 1 illustrates this: While the surface patch on the left is locally homeomorphic to Euclidean space, the one in the middle is a Tjunction, and the one on the right is a touching point. Twodimensional domains with such features appear in applications like the simulation of closedcell foams (Nammi et al., 2010), or networks of fractures in rock mechanics (McClure and Horne, 2013; McClure et al., 2015). Of considerable importance are also onedimensional networks embedded into a two or threedimensional Euclidean space. They appear in models of traffic networks (Garavello and Piccoli, 2006), supply chains (D’Apice et al., 2010), but also in simulations of biological systems like root networks (Dunbabin et al., 2013)
(Lang et al., 2011), or blood vessel networks (Cattaneo and Zunino, 2014a). An overview over flow problems on networks is given in (Bressan et al., 2014). Figure 2 shows two example domains for network problems.Various discretization methods have been proposed for surface and network equations (Dziuk and Elliott, 2013; Olshanskii et al., 2009). Explicit discretizations, which are the focus of this work, use a grid of the same dimension as the domain. For manifold surface grids it is reasonably simple to generalize grid data structures to such a setting. The main hurdles are admitting that the number of coordinates of a vertex can be different from the effective grid dimension, and making sure that grids are not required to have a boundary. Several standard simulation codes support such surface grids. We mention Alberta (Schmidt and Siebert, 2005) (standalone and as part of Dune), AMDiS (Vey and Voigt, 2007), and FEniCS (Rognes et al., 2013). Moreover, the GeometryGrid Dune meta grid allows to embed any Dune grid into a Euclidean space of higher dimension.
Grid data structures for nonmanifold grids are more challenging. To handle Tjunctions, for example, the data structure must cope with the fact that element facets (i.e., edges in a 2d grid) may have more than two neighbors (Figure 1). While this is not very difficult to implement, it requires the introduction of additional data fields and logic, which is not used when the grid happens to be a manifold. Since the latter case is predominant, such additional features mean space and runtime overhead for most users. Therefore, standard grid data structures do not allow for nonmanifold topologies. Numerical examples of processes on nonmanifold topologies in the literature are typically done using ad hoc implementations of the necessary grid data structures.
Unfortunately, using ad hoc implementations wastes a lot of human resources. While the domain may be nonstandard, many of the equations on network grids differ little from their Euclidean counterparts. However, existing code like finite element assemblers for diffusion and transport processes cannot be reused on ad hoc grid data structures. They require detailed knowledge of the grid, and it is a challenge to port an existing assembler to a new grid data structure. Also, since the data structure is ad hoc, it is difficult to reuse it in other contexts. In particular, it is difficult to share it with other groups working in the same field.
The Dune software system^{1}^{1}1www.duneproject.org (Bastian et al., 2008b, a) has found an elegant solution for both problems. Dune is a set of opensource C++ libraries dedicated to various aspects of finite element and finite volume methods. Its grid component, implemented in the dunegrid module, specifies an abstract interface for computational grids. The specification mandates what a data structure should be able to do to qualify as a Dune
grid, and how this functionality should be accessible from C++ code. Examples of such functionality include being able to iterate over the elements and vertices, and getting the maps from the reference element to the grid elements. Those parts of a numerical simulation code that use the grid, such as the matrix assemblers or error estimators, are written to use only the abstract grid interface. Grid data structures implementing this interface are then completely decoupled from the algorithms that use them.
This decoupling has various interesting consequences. Using the Dune grid interface, it is easy to swap a given grid implementation for another one. Indeed, in typical Dune applications the C++ grid type is set once in the code, and handed around as a template parameter. Changing the initial typedef and recompiling the code is usually sufficient to switch to an alternative grid data structure. This possibility to easily switch between grid implementations allows to provide tailormade grid data structures for special simulation needs. For example, dunegrid itself provides YaspGrid, the implementation of a structured grid with very little runtime and space overhead. In contrast, UGGrid implements a very flexible unstructured grid with nonconforming and red–green refinement.
In this paper we present FoamGrid, a new implementation of the Dune grid interface that is dedicated to surface and network grids. Its main features are:

FoamGrid implements one and twodimensional simplex grids embedded in a physical Euclidean space of arbitrary dimension . Hence, it targets geometric and surface PDEs.

A FoamGrid can be adaptively refined. In the standard setup, refining a triangle results in four coplanar triangles. Additionally, FoamGrid elements can be parameterized, i.e., they can be given a map that describes a nonlinear embedding of the element into . As the element gets refined more and more, its shape approaches the one described by the parameterization (Figure 4).

Finally, the domain of a FoamGrid can grow and shrink at runtime. Elements can be added and removed even when there is no coarser “father” element, without invalidating the grid data structure. Data can be transferred during this process. This allows to elegantly simulate network growth and remodeling processes.
FoamGrid is available in a Dune module dunefoamgrid,^{2}^{2}2http://users.duneproject.org/projects/dunefoamgrid and is installed just like any other Dune module. dunefoamgrid is free software, available under the either the LGPLv3+, or the GPLv2 with a linking exception clause.
Using Dune and FoamGrid for simulations of network and surface PDE problems has a number of important advantages. First of all, you do not have to implement the grid data structure yourself. While not overly difficult, quite a bit of thought has gone into FoamGrid, which would need work to be replicated. Then, since FoamGrid implements the Dune grid interface, large amounts of existing application and infrastructure code can be used directly with FoamGrid. This includes things like finite element spaces and assemblers, error estimators, and grid file readers and writers. This advantage is demonstrated in particular by the numerical example in Section 4.1, which required no additional coding at all to extend an existing planar code to a network setting.
As a further advantage, once a user starts employing FoamGrid and the Dune grid interface, he immediately has the power of all other Dune grid implementations at his disposal. All of these are easily usable together with FoamGrid. Hence, e.g., in simulations that couple network grids with background grids, several implementations of such background grids are readily available.^{3}^{3}3And the Dune extension module dunegridglue offers a convenient way to do the coupling (www.duneproject.org/modules/dunegridglue). Since the Dune grid interface is a wellestablished standard, it is easy to learn how to use FoamGrid. If a user already knows Dune, there is little additional knowledge needed. Since FoamGrid is open source, sharing code based upon it is particularly easy.
While FoamGrid has many interesting features, there is also a number of things it does not currently support. For example, elements can only be simplices, and they must be one or twodimensional. (The authors could not think of a use case for a higherdimensional network grid. Write us if you know one.) Adaptive grid refinement is currently nonconforming, which leads to hanging nodes, and, possibly, to holes in the surface (Figure 5). Finally, the current FoamGrid implementation is purely sequential, and FoamGrid objects cannot be distributed across several processes. However, the development of FoamGrid is ongoing, and these features may appear in later releases.
The present article is structured as follows. Chapter 2 briefly explains how to use FoamGrid. Everything mentioned there is codified in the Dune grid interface specification, and hence you can also read this chapter as an introduction to the use of the Dune grid interface. Chapter 3 then explains the special features that FoamGrid offers. In particular, these are support for Tjunctions, adaptive refinement with element parameterizations, and the ability to “grow”. Finally, in Chapter 4 we give three numerical examples showcasing the different features of FoamGrid. The first shows unsaturated flow through a twodimensional fracture network using finite elements. The second example shows adaptive, locally massconservative transport of a therapeutic agent in a microvascular network using finite volumes. The third one models the growth of plant root networks.
2 FoamGrid and the Dune grid interface
In this chapter, we start by describing the programmer interface of FoamGrid. FoamGrid implements the Dune grid interface, hence in many central aspects, it can be used just like any other Dune grid. This chapter focuses on these aspects. You can therefore also read it as a brief review of the Dune grid interface. For more details, you may want to consult the Dune online documentation and (Bastian et al., 2008b, a).
The central class of the FoamGrid grid implementation is
available from the header dune/foamgrid/foamgrid.hh. This class implements a hierarchical grid as defined in (Bastian et al., 2008b, Def. 13), i.e., a coarse (or macro) grid, and element refinement trees rooted in each of its elements. The first template parameter dim is the grid dimension , which must be either 1 or 2. The second template parameter dimworld is the dimension of the Euclidean embedding space. It must be equal to or greater than the grid dimension, but can otherwise be arbitrary. For the rest of this article we use dim and dimworld in code examples, and and in text to denote the dimensions of the grid and the physical space, respectively.
To construct FoamGrid objects, the FoamGrid class implements the entire Dune grid interface for the setup of unstructured grids. In particular, the class GridFactory is implemented for FoamGrid. Thus, all filereading methods based on this interface are available. For example, files in the gmsh format (Geuzaine and Remacle, 2015) can be read by using the line
This will read the file named “filename.msh”, and set up a new FoamGrid<2,3> object with it. The new grid object is returned in a shared pointer called grid. Note that vertex coordinates in gmsh files always have three components, so reading a gmsh file into a FoamGrid<1,2> object will discard the third entry. As a special feature, gmsh files can contain elements with polynomial geometries of order up to five. While FoamGrid element geometries are always affine, FoamGrid can use the higherorder geometries during mesh refinement (Section 3.2).
Writing FoamGrid objects to disk is equally straightforward. All Dune file writing codes rely on the grid interface only, and can therefore be used with FoamGrid. For example, writing the object pointed to by the grid shared pointer into a VTK file called my_filename.vtu can be achieved by including the header dune/grid/io/file/vtk.hh from the dunegrid module, and writing
The resulting file can be visualized, e.g., with the ParaView software.
2.1 Elements and geometries
In Dune, finite element assembly typically takes place on the leaf elements of the refinement trees. These elements form the leaf grid view, which encapsulates the notion of a textbook nonhierarchical finite element grid. From a FoamGrid, the leaf grid view can be obtained in the usual way, i.e.,
which already has been used in the VTK example above.
Access to grid elements and vertices is provided by the grid view. Elements can be iterated over using begin/enditerators
where the number specifies that the loop is to be over the grid elements (it is the codimension of the grid elements with respect to the grid). Using the C++11 rangebasedfor syntax, the same loop can be written more concisely
Similarly, a loop over all vertices of the grid is written as
In this code, the number dim (i.e., the grid dimension), specifies that the loop is to be over the grid vertices, because vertices are zerodimensional and hence have codimension dim in a dimdimensional grid. Alternatively, one can write
The objects vertex and element (or *it in the iterator loops) are instances of what in Dune terminology are called entities. Entities are implemented by the Entity interface class in dunegrid. They provide topological information about the grid elements and vertices, like links to the corners vertices of an element, and to the father and descendant elements in the refinement tree. In all these aspects, FoamGrid objects behave just like any other Dune grids. Furthermore, each element entity provides a Geometry object, which represents the affine map from the reference element to the grid element . This map provides the geometrical information needed to assemble finite element and finite volume systems, like evaluation of and its inverse , the inverse transposed Jacobian matrix , and the functional determinant . Note that since for surface grids the world dimension is larger than the dimension of the reference element , the image of under is a set of measure zero in . In finiteprecision arithmetic the argument of the method implementing will typically not be in the domain of . The Geometry implementation of FoamGrid therefore extends to the entire space . Given a point not necessarily in , FoamGrid computes a point in the plane spanned by such that is minimized. The affine function is tacitly extended from to its entire affine hull here.
To compute element fluxes, elements in a Dune grid provide a set of socalled intersections, which relate elements to their neighbors. This mechanism is very flexible, and in particular handles general nonconforming situations very well. Nevertheless, using grids for network domains stretches the bounds of the current intersection concept, and some generalization is needed. We have therefore dedicated a separate chapter to FoamGrid intersections (Chapter 3.1).
2.2 Attaching data to grids
Data is attached to FoamGrid objects in the same way as to other Dune grids. Each grid view object can provide a corresponding IndexSet object
This index set provides an integer number for each entity (i.e., vertex, edge, or element) of the grid view. For each dimension and reference element type, these numbers are consecutive and start at zero. They can hence be used to address randomaccess containers holding the simulation data. This approach is convenient, flexible, and efficient.
Data stored in arrays is lost if the grid changes, either by refinement (Section 2.3) or by grid growth (Section 3.3). To preserve data across grid modifications, FoamGrid, just like any other Dune grid, additionally provides a set of persistent numbers. These are obtained by an IdSet object
(remember that in our initial example the variable grid was a shared pointer to a FoamGrid). Persistent numbers are neither consecutive nor restricted to start at zero, but they can be used to access search trees or hash maps. Before modifying the grid, all simulation data must be copied into such data structures, and copied back to arrays after the modification is completed. While such copying is costly, its runtime is usually negligible compared to the cost of the actual grid modification.
2.3 Adaptive refinement
FoamGrid supports red refinement (nonconforming refinement) of simplices, where each triangle is split into four congruent smaller triangles. If a twodimensional FoamGrid is refined locally, then hanging nodes appear in the grid. Depending on the discretization used, this may or may not be a problem. True red–green refinement, which avoids the hanging nodes, may appear in later versions of FoamGrid.
Adaptive grid refinement in FoamGrid is controlled via the standard Dune grid interface. In a first step the method
is used to mark an element element for refinement (refCount ) or coarsening (refCount ). The method returns true if the element was successfully marked. The mark of an element element can be obtained with the method
which returns 1 for elements marked for refinement, 1 for elements marked for coarsening, and for unmarked elements.
The grid is then modified in a second step with the methods
Between preAdapt() and adapt() it is possible to check the following flag
Similarly, between adapt() and postAdapt() one can check the flag
Both are useful to manage the transfer of data associated with the grid. As FoamGrid itself does not store any associated data, the data transfer from the old grid to the adapted grid is managed by the user. An example featuring adaptive refinement and coarsening is presented in Section 4.2.
3 Additional features
The previous chapter has described those aspects where FoamGrid behaves just like any other Dune grid. However, FoamGrid also has a few features that set it apart from most other Dune grids. The present chapter is dedicated to those features.
3.1 Handling intersections in a nonmanifold grid
The defining feature of FoamGrid is its capability to handle network grids, i.e., grids with a nonmanifold topology (Figures 1 and 2). In particular, more than two elements can meet in a common vertex in a onedimensional grid, and more than two triangles can meet in a common edge in a twodimensional grid.
Information about how a given element relates to its neighbors is essential for finite volume and DG methods, or more generally all methods involving element boundary fluxes. For this, the Dune grid interface provides the notion of intersections. In Dune terminology, an intersection is the settheoretic intersection between the closures of two neighboring elements. Only intersections that have positive –dimensional measure qualify as intersections in the Dune sense, and those are equipped with a coordinate system, i.e., a map from a –dimensional reference element to the intersection (Figure 3). Note that if the grid is conforming, then intersections will be geometrically the same as shared element faces. However, in general nonconforming grids, intersections are only subsets of element faces.
Intersections provide all information needed to compute element boundary fluxes. In particular, for each point on an intersection, you can get the vector
that is tangent to the element at that point, and normal to the element boundary. Also, you get the embeddings of the intersection into the two elements, in form of maps from the intersection reference element to the element reference elements (Figure 3). Those maps pick up the same ideas used to model element geometries in Section 2.1. They are called geometryininside and geometryinoutside, respectively.In C++ code, intersections are objects of type Intersection. For any given element (which is then called the inside element), all its intersections with neighboring elements can be accessed by traversing them with a dedicated iterator. Using rangebased for syntax, a loop over all intersections of the element called element is written as
See the dunegrid class documentation for the details on the user interface of the Intersection class.
Unfortunately, the Dune intersection mechanism, as flexible as it is, is not flexible enough for network grids. The original idea was that while elements may have more than one neighbor across a given facet, there is (even in nonconforming situations) at most one neighbor at any given point on that facet. This reflects the assumption that computational domains are expected to be topological manifolds. At the time of writing, this assumption is still reflected in the grid interface. In particular, the method
(a public member of the Intersection class), informs whether there is a neighbor across a given intersection. However, this information is insufficient in network grids, where even in a conforming grid there may be an arbitrary number of neighbors meeting at a common intersection. Finite volume and DG methods need access to this group of neighbors, to know how to distribute the flux across this intersection.
For this reason the intersection concept and programmer interface is being revisited, and is likely to undergo changes in the future to better support network grids. Changes to the Dune interface can only be made in a democratic process, and it is therefore unclear at what point in time such a change will happen.
Two different approaches have been proposed to the Dune grid development community. We describe them both briefly here. The full proposal text is available at www.duneproject.org/modules/dunefoamgrid.
The first approach tries to be as minimally invasive as possible. In particular, it retains the notion of an intersection as an object that relates two elements. The extension consists of two semantic rules, and a change to the neighbor method. The first of the semantic rules makes sure that the “number of neighbors” across a given intersection is welldefined. Remember that the geometryininside is, roughly speaking, the intersection interpreted as a subset of the element .
Rule 1
For any two intersections of a given element , the geometriesininside are either disjoint or identical.
With the number of neighbors properly defined, the neighbor method is generalized to return this number instead of a yes/no answer:
Note that this change is fully backwardcompatible, as intersections in a nonnetwork grid will return either 1 or here, which casts to the values true or false as used previously.
The second semantic rule provides a way to find all intersections that share a common geometryininside. No additional interface method is added for this. Rather, it is guaranteed that all such intersections appear consecutively when traversing the intersections with the intersection iterator.
Rule 2
If more than one neighbor is reachable over a given geometryininside, then all intersections for this geometryininside shall be traversed consecutively by the intersection iterator.
With this rule, groups of neighbors can be identified and used in flux computations.
The second approach is more radical, because it changes the idea of an intersection itself. Intersections cease to be objects that relate pairs of elements. Rather, they now become objects that relate groups of elements. As a consequence, each intersection still has only one geometryininside. However, for each intersection there is now more than one outside element, each with corresponding geometryinoutside and indexinoutside.
To access this information, more interface methods need to be changed. First of all, the neighbor method needs to be changed as proposed above. Secondly the methods
of the Intersection interface class need to be replaced by
respectively. Rather than returning the unique outside element or its geometry or index, the methods must now return the corresponding quantity for the ith outside element.
These changes are again fully backwardcompatible. In grids without multiple intersections, at most the th outside element will be available. The default parameter ensures that this intersection will be returned when the method is called without argument.
As an advantage, this proposal retains the rule that the geometriesininside must form a disjoint partition of the element boundary (modulo zerosets). Also, it is easier to attach data to such intersections, which is a feature that has been requested various times in the past. On the downside, to iterate over all neighbors of an element, two nested loops are needed, instead of only one. This will make some code a bit longer, and more difficult to read.
Both proposals are currently under discussion. However, even with the current status quo, applications involving fluxes can be written for network domains. Onedimensional domains are straightforward as there the intersections are a fortiori conforming (see Sections 4.2 and 4.3). Twodimensional networks need more trickery, but can also be made to work. Once either of the proposed interface extensions has been officially accepted, implementations of such methods will be much simpler.
3.2 Element parametrizations
In a standard nonconforming red refinement algorithm, new vertices are placed at the edge midpoints of refined elements. That way, while the grid gets finer and finer, the geometry of the grid remains identical to the geometry of the coarsest grid. To also allow improvements to the geometry approximation, FoamGrid can use element parametrizations. Each coarse grid element with reference element of a FoamGrid can be given a map
which describes an embedding of into physical space . This does not influence the grid itself — elements of a FoamGrid are always affine. However, when refining the grid, new elements are not inserted at edge midpoints. Rather, for each new vertex which would appear at an edge midpoint in standard refinement, the corresponding local position in its coarsest ancestor element is determined using the refinement tree. This position is then used as an argument for , and the result is used as the position of the new vertex. That way, the grid approaches the shape described by the parametrization functions more and more as the grid gets refined (Figure 4).
In dunegrid, element parametrizations are implemented as small C++ objects. These must inherit from the abstract base class
declared in dune/common/function.hh.^{4}^{4}4For various reasons there is discontent with using the VirtualFunction class to implement element parametrizations. Therefore, it is expected to be replaced by something else eventually. However, at the time of writing, there is no specific proposal for such a replacement. In any case, while the details of the implementation are under discussion, the general idea of element parametrizations is not disputed, and can be expected to remain the way it is described here. One object of this type needs to be created for each element of the coarse grid. The base class has a single pure virtual method
which implements the evaluation of . The first argument x is a position in local coordinates of the corresponding coarse grid triangle. The result is returned in the second argument y.
Element parametrizations are handed to the GridFactory during grid construction. Normally, grid elements are entered using the method
of the GridFactory class. For parametrized elements, there is the alternative method
This inserts the element with vertex numbers given in the array vertices and an element parametrization given by the object elementParametrization into the grid. The GmshReader uses this method for some of the higherorder elements that can appear in gmsh grid files.
To see the effect of element parametrizations, Figure 4 shows an example where the coarsest grid consists of only two triangles covering the domain . As a parametrization we use the global function
(1) 
and each element parametrization first maps its local coordinates to , and then applies the global function . Hence, upon refinement, the grid will approach the graph of the function (1). The code for this example is provided in the dunefoamgrid module itself, in the file dunefoamgrid/examples/parametrizedrefinement.cc.
There is one pitfall when using element parametrizations together with nonconforming adaptive refinement for twodimensional grids. If a hanging node appears in the course of refinement, the position of this node is determined by the parametrization functions of the corresponding element, which is typically not the midpoint of the adjacent longer edge. As a consequence, the grid will have holes wherever different refinement levels meet (Figure 5). While this may seem surprising at first sight, it is nevertheless a logical consequence of the hanging nodes refinement. The continuous domain is approximated by discontinuous grids, in direct correspondence to how discontinuous FE functions may be used to approximate continuous PDE solutions.
The nonconforming geometry approximation shows another shortcoming of the current Dune intersection concept. There, intersections are defined as settheoretic intersections of (closures of) neighboring elements. However, in the geometrically nonconforming situation, neighboring elements do not actually intersect, and one can only speak of logical intersections. The practical consequence of this is that intersections do not have a single welldefined shape in anymore. Rather, there are now two of them. This possibility to have two global element shapes is not currently reflected by the Dune grid interface. In the current implementation of FoamGrid, the method Intersection::Geometry will always return the global shape of the intersection as seen from the inside element.
3.3 Grid growth
A certain number of network problems is posed on domains that grow and/or shrink in the course of the simulation. Examples are fracture growth processes and simulations of bone trabeculae remodeling. To support such simulations, FoamGrid objects are allowed to grow and shrink, i.e., elements can be added and removed from the grid at runtime. Finite element and finite volume data is kept during such grid changes, using an approach much like the one used for grid adaptivity. Of all Dune grids, FoamGrid is currently the only one to support this feature. The programmer interface combines ideas from the GridFactory class to insert new vertices and elements, with the adaptivity interface to allow to keep data across steps of grid growth.
Growth and shrinkage of grids is a twostep process. First, new elements and vertices are handed to the FoamGrid object. These are not inserted directly; rather, they are queued for eventual insertion. In addition, individual elements can be marked for removal. Once all desired elements and removal marks are known to the grid, the actual grid modification takes place in a second step.
Queuing elements for insertion and removal is controlled by three methods. The first,
queues a new vertex with coordinates x for insertion. The return value of the method is an index that can be used to refer to this vertex when inserting new elements. The index remains fixed until all queued elements are actually inserted in the grid (by the grow method), but may change during the execution of that method. Inserting elements is done using
which mimics the corresponding method from the GridFactory class. The argument type has to be a simplex type, because (currently) FoamGrid supports only simplex elements. The array vertices must contain the indices of the vertices of the new element to be inserted. These can be either indices of existing vertices, or new indices obtained as the return values of the insertVertex method.
Analogously a new element with a parametrization can be inserted by calling
Finally, the method
marks the given element for removal.
Once all desired elements are queued for insertion or removal, the actual grid modification takes place in a second step. The grid is modified using the method
While element removal is guaranteed, queuing elements does not assure that the element will be inserted. New elements are restricted by the fact that Dune grids are hierarchic objects. The vertices given by the user to form an element are always leaf vertices but may be contained in different hierarchic levels. However, elements can only be constituted by vertices of the same level. Therefore, new elements in FoamGrid are always inserted on the lowest possible level substituting the given vertex by its hierarchic descendants or ancestors. Note that it is not generally guaranteed that relatives of the given vertices on the same level can be found. In that case, the element will not be inserted. The method grow will return true if it was possible to insert at least one element.
After the call to grow, it is possible to check whether a given element has been created by the last call to the grow method:
which is a method of the interface class Entity<0>, i.e. elements. Observe that this is the same method that returns whether an element has been created by grid refinement. Hence its semantics depends on whether it is queried after a call to grow or after a call to adapt. Using this method is helpful, e.g., when setting initial values and/or boundary conditions for newly created elements.
The growth is completed with the call
which removes all isNew markers.
Summing up, growing or shrinking the grid while keeping grid data consists of the following steps. Note the relationship to grid adaptivity with data transfer.

Mark elements for removal; queue new vertices and elements for insertion.

Transfer all simulation data attached to the grid into an associative container indexed by the entity ids described in Section 2.2.

Call grow().

Resize data array; copy data from the associative container into the array.

Set initial data at newly created elements and vertices and boundary conditions at newly formed boundaries. Note that element removal always creates new boundaries.

Finalize by calling postGrow().
While grid growth itself is straightforward, it is difficult to use in combination with adaptive refinement. Grid refinement in Dune leads to hierarchical grids, which are forests of refinement trees. Not every element of such a forest can be added or removed without violating certain consistency conditions. In onedimensional grids, there are relatively few problems, and grid growth and refinement can be used together to good effect. For twodimensional grids we have tried to be as general as possible, but there are limits.
There are obstacles both to the removal of elements and to the insertion of new ones, if grid refinement is involved. First of all, only leaf grid elements can be removed. There is no conceptual problem with element removal if the element has no father. If it does have a father, however, removing only a subset of its sons is a violation of the Dune grid interface specification, which mandates that the sons of an element must (logically) cover their father (Bastian et al., 2008b, Def. 11.1). Deliberately allowing this violation nevertheless appears to be the only viable solution, as all other possibilities amount to only allowing the removal of large groups of elements, which is rarely desired.^{5}^{5}5Another Dune module which violates this assumption to good effect is dunesubgrid, see (Gräser and Sander, 2009).
Inserting new elements into a refinement forest of elements is even more problematic, because the new element needs to be assigned a level number in the hierarchy. Ideally, this level would always be zero, because if the element had a larger level number it would be expected to have a father. FoamGrid does violate this assumption if necessary, which does not lead to problems in practice, unless multigridtype algorithms are used. On the other hand, the element level must be the same as the levels of all of its vertices. If the vertices are new, their levels can be freely chosen. However, grid growth almost always involves vertices that already exist in the grid. When trying to insert elements with vertices having different level numbers, FoamGrid currently tries to replace existing vertices by their father vertices, to obtain a set of vertices on the same level (which then determines the element level).
4 Numerical examples
We close the article with three numerical examples. These show how seemingly challenging algorithms can be implemented with ease using the FoamGrid grid manager.
4.1 Unsaturated Darcy flow in a discrete fracture network
For our first example we simulate unsaturated Darcy flow through a network of twodimensional fractures embedded into . We only consider the flow in the network itself, but FoamGrid can be easily coupled to higherdimensional background grids using the dunegridglue module (Bastian et al., 2010; Engwer and Müthing, to appear).
Let the computational domain be the union of a finite number of closed, bounded hypersurfaces in . We use the Richards equation to model the flow in (Bear, 1988). That is, we suppose that the state of the system in a time interval can be described by single scalar pressure field
From this pressure field, the fluid saturation and permeability can be computed using the Brooks–Corey and Burdine parameter functions
where , , , and are scalar parameters. The saturation satisfies the Richards equation
(2) 
For simplicity we suppose that the flow is purely driven by the boundary conditions, and omit the gravity term .
Equation (2) is a quasilinear equation in the pressure . In (Alt and Luckhaus, 1983) (see also (Berninger et al., 2011)) it was shown how the Kirchhoff transformation can be used to transform it to a semilinear equation for a generalized pressure
We discretize this equation in time using an implicit Euler method, and in space using firstorder Lagrangian finite elements. The resulting weak discrete spatial problem can be written as a minimization problem for a strictly convex functional. At each time step, we determine the minimizer of this functional by a monotone multigrid method (Berninger et al., 2011).
For the implementation we used the code used for the numerical examples in (Berninger et al., 2011). Since vertexbased finite elements were used for the discretization, no changes to the numerical algorithm were needed to also apply it to network grids. Originally, the code used the Dune libraries and the UGGrid grid manager for unstructured grids. Even though the code for (Berninger et al., 2011) was not written with network flow problems in mind at all, it could nevertheless be reused as is, after only a handful of trivial bugfixes. The only changes necessary were replacing the line
by
and adjusting the boundary data specification. This once more proves the point that using the Dune grid interface gives great flexibility, and allows code to do more things than planned by the original authors.
We present an example simulation using an artificially created network grid. The grid, which can be seen in Figure 6, was created by Patrick Laug using the fracture network grid generator described in (Borouchaki et al., 2000). The grid spans the volume (length and pressure are given in meters). We assume the network to be filled with a sandlike material with parameters , , , bubbling pressure m, and .
We assume the network to be initially devoid of water, setting m. Then, water is injected in a unit circle centered at the point on the boundary with a constant pressure of m; noflow boundary conditions are imposed at the remaining boundary. Figure 7 shows several steps in the evolution of the physical pressure . As expected, one can see the fluid entering, and slowly filling almost the entire network. At the end, a steady state is reached, and the pressure is constant in each connected component of the domain.
4.2 Transport of a therapeutic agent in the microvasculature
The next example demonstrates local grid adaptivity, and the treatment of network bifurcations. We model the propagation of a therapeutic agent for cancer therapy in a network of small blood vessels. To reduce the computational effort, the network of threedimensional vessels is reduced to a onedimensional network embedded in a threedimensional tissue domain .
We first discuss the onedimensional model for flow in a blood vessel segment, disregarding any bifurcations. Starting from a threedimensional blood vessel domain , we describe blood in the microvasculature as an incompressible Newtonian fluid with viscosity and density governed by the Stokes equation. Further, we assume an axially symmetric blood vessel segment with constant radius , crosssection area , the parametrized tangent on the vessel centerline , , , and a rigid vessel wall. Then, with the assumption of constant pressure in a crosssection and negligible radial velocities the Stokes equations can be integrated over a vessel segment yielding the onedimensional equation
(3) 
The parameter shapes a powertype axial velocity profile with mean velocity , yielding a quadratic velocity profile for and flatter profiles for . For a detailed derivation in a more general setting, see the reduction of the Navier–Stokes equations to onedimensional equations in (Quarteroni and Formaggia, 2003). Equation (3) is a simplified stationary version of the onedimensional blood flow equations described by Quarteroni and Formaggia (2003), neglecting vessel wall displacement and intertial forces.
Small blood vessels are exchanging mass with the embedding tissue through the vessel wall. The fluid exchange can be modeled by Starling’s law, and results in an additional source term. With this modification, (3) becomes
(4) 
where is the empirical filtration coefficient dependent on, e.g., the intrinsic permeability and thickness of the membrane, and the viscosity of the interstitial fluid. The source term further depends on the pressure in the surrounding tissue . For the sake of simplicity, this tissue pressure is subsequently assumed constant. Equation (4) was also used by Cattaneo and Zunino (2014a) in a finite element setting to model coupled vesseltissue flow processes in a tumor tissue.
The transport of a therapeutic agent is modeled by an advection–diffusion equation using the velocity field calculated by equation (4). Similar to the reduction of the Stokes equations we can reduce the threedimensional advection–diffusion equation by integration over a vessel segment assuming a constant concentration on a given crosssection with area (D’Angelo, 2007; Cattaneo and Zunino, 2014b). The transport over the vessel wall can be described by the Kedem–Katchalsky equation (Kedem and Katchalsky, 1958), yielding
(5) 
The last term in (5), accounting for transport across the vessel wall, consists of an advective and a diffusive part where advection is reduced by the reflection coefficient for larger molecules. Again, we assume the mean tissue concentration to be constant.
To model a network of such segments we split the blood vessel network at junctions into pieces yielding a set of vessel segments each governed by equations (4) and (5). At each junction we require continuity of pressure
Together with these coupling conditions, and boundary conditions on , Equation (4) has a unique solution, which is the stationary pressure field .
For the transport at the junctions, we require continuity of concentration
and, for , the volume flux leaving segment at the junction, we require mass conservation
We discretize equations (4) and (5) in space with a standard finite volume scheme and piecewise linear onedimensional grid elements. This demands balancing fluxes over the edges of the elements. Assume that we have calculated all element transmissibilities
(6) 
Then, the volume flux from element to a neighboring element can be calculated as
(7) 
One can see that the transmissibility at branching points is dependent on the transmissibility of all neighboring elements. Assuming that all element transmissibilities reside in an array transmissibility, the calculation of the twopoint transmissibilities could look as follows using FoamGrid. Note that only a single loop over the grid entities is necessary to calculate the transmissibilities.
This code works well using the grid interface of dunegrid2.4, even though that interface does not have any provisions for network grids at all (Section 3.1). This is because the grid is onedimensional. In this case, each intersection always covers entire element facets (i.e., vertices). Therefore, the indexInInside method can be used to group intersections.
To reduce numerical diffusion induced by the implicit Euler scheme, the grid can be refined adaptively around the concentration front. Initially, the grid is refined around inflow boundaries, and during the simulation the grid is adapted by a local gradientbased concentration indicator as described in (Wolff, 2013). This indicator marks an element for refinement if
and for coarsening if
where denotes the concentration difference between element and a neighboring element . The parameters , () are problem dependent. We achieved robust adaptation behavior in our example for and . The indicator is evaluated before every time step and marks elements for refinement or coarsening using the mark method, see Section 2.3. The adaption of the grid is then handled by FoamGrid. Between the preAdapt and postAdapt steps, we have to transfer data, i.e., primary variables and spatial parameters from the old grid to the new adapted grid. As can be seen in Figure 9, the refinement scheme works well even around network bifurcations.
We simulate a network of capillaries in rat brain tissue scanned by Motti et al. (1986) and reconstructed as segment network with threedimensional geometrical information by Secomb et al. (2000). The network data comprises threedimensional location of vessel segments, inflow and outflow boundary markers, vessel segment radii, and velocity estimates for each vessel segment. The domain has a bounding box of m. The given vessel radii vary between m and m. The estimated velocities are in a range of to . We choose the blood viscosity , the filtration coefficient , the effective diffusion coefficient over the vessel wall , and the diffusion coefficient of the transported agent trail , calculated with the Stokes–Einstein radius. We assign the following boundary conditions for the flow problem (4)
The velocities for the inflow segments are set to the estimates provided by Secomb et al. (2000) along with the grid geometry. We simulate the arrival of a therapeutic agent by a Dirichlet boundary condition for equation (5) on a subset of , namely,
Specifically, we enforce a mole fraction of at one of the inflow vessel segments with a Dirichlet boundary condition. At outflow boundaries, we neglect diffusive fluxes. Note that a full upwind scheme is employed for the concentration, so no further boundary condition for the advective fluxes is necessary at outflow boundaries.
Figure 8 shows the resulting stationary pressure field and the initial element size. In Figure 9, one can see the resulting mole fraction at various times. Note how the local grid refinement follows the steepest gradients, i.e., the transport front. The resulting mole fractions vary from segment to segment due to different radii. This also automatically ensures a finer grid around vessel bifurcations. The onedimensional elements are depicted as threedimensional tubes scaled with their respective radius.
4.3 Root water uptake and root growth at plantscale
In environmental and agricultural research fields, models describing root architectures are used to investigate water uptake and root growth behavior of plants (Dunbabin et al., 2013). In our final example, a onedimensional network embedded into is used to describe such a plant root architecture. We simulate water flow through the root network and root growth.
Plant roots can be described by a treelike network of pipes which consist of xylem tubes (Tyree and Zimmermann, 2002). We follow the cohesion–tension theory (Tyree, 1997), where the water flow through the root system is governed by the pressure gradient caused by the transpiration rate of the plant above the soil. We assume only vertical flow and no gravity. This leads to a Darcy’s law analogy (Doussan et al., 1998)
where is the water flux in the xylem tubes, is the xylem water pressure, and is the axial conductance of one root segment.
Water can enter the roots at any point on the xylem tube surfaces, which leads to a volume source term for the onedimensional network model. For simplicity, we model a single membrane only for the entire pathway of water from soil into the roots. With this assumption, and neglecting osmotic processes, radial water flow into one root segment is defined as
where is the radial conductivity, i.e., the conductivity of series of tissues from root surface to the xylem. The number is the soil–root interface area. The water pressure at the soil–root interface must be provided. One option is to couple the root system to a Richards equation based soil water flow simulation (Javaux et al., 2008), but for simplicity we simply take as a known value.
The continuity equation leads to
(8) 
with a solutiondependent source
where is the only unknown variable. This modeling approach neglects the influence of solutes on water flow, as well as the capacitive effect of the roots, because the amount of water stored in roots is generally small compared to transpiration requirements.
Equation (8) is discretized in space with standard cellcentered finite volumes, piecewise linear onedimensional grid elements and implemented using the external Dune discretization module DuMu (Flemisch et al., 2011), and the FoamGrid grid manager. Flux calculations over edges and branching points of the root network are implemented just as in our previous example.
So far, we have assumed a root network that does not alter its geometry over time. Several algorithms were developed to describe root growth (e.g., (Pagès et al., 2004; Somma et al., 1998; Leitner et al., 2010)). These models define root growth either by a fractal description, or more generically depending on plant specific parameters (root elongation, growth direction, branching density) and surrounding soil properties (soil moisture, soil strength, temperature, nutrients). For simplicity, we model root growth here as an (almost) completely random process. New root branches occur at random time steps and with a small gravity effect only, which makes the roots tend to point downwards. Existing branches grow at the branch tip at random times and without changing directions.
Our example simulation starts with a simple root grid which consist of one vertical root branch discretized with eight elements (root segments) and no lateral branches (Figure 10). We choose radial and axial conductivity values from (Doussan et al., 1998). The surrounding relative soil pressure is set to Pa, and the Dirichlet boundary value at the root collar is set to Pa. Parameters and boundary conditions do not change with time.
In every time step, new root elements are created and either added to an existing lateral branch or to the main branch. An elementbased indicator based on simulation time and random factors decides whether a new branch is added at one of the element’s vertices. The Indicator class also computes the coordinates of the newly inserted point that is the second vertex of the new element. The coordinates depend again on simulation time, random factors, and the branch orientation in . Our growth step calculation in DuMu using FoamGrid looks as follows.
The vars container contains all primary variables and spatial parameters defined on the root segments. The Indicator class is a template parameter that can be easily exchanged, to allow different root growth algorithms. In Step (2), the new elements are inserted. New root segments must be connected to the old grid. The implementation of the insertElements method could looks as follows:
The primary variables and spatial parameters of our physical problem are stored in Step (3), before the actual growth step (4). We update the sizes of the variables and parameter vectors since the total number of degrees of freedom has changed due to the growth step (5). In addition, values for the new elements have to be computed. In our case, new root segments inherit the primary variables and the spatial parameters from its preceding neighbor element. Root tips, in particular new root tips, are always assigned Neumann noflow boundary conditions. At the end, the
postGrow method is called, which deletes the isNew markers (Step (6)).Figure 10 shows the root network and the pressure distribution inside the roots for several time steps. The total root water uptake (transpiration demand) of the plant, defined by the Dirichlet boundary condition, does not change during the growing period. Thus, the pressure inside the plant changes since the water uptake of the plant is distributed to more and more roots segments.
References
 Alt and Luckhaus [1983] H. Alt and S. Luckhaus. Quasilinear elliptic–parabolic differential equations. Math. Z., 183:311–341, 1983.
 Bastian et al. [2008a] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, R. Kornhuber, M. Ohlberger, and O. Sander. A Generic Grid Interface for Parallel and Adaptive Scientific Computing. Part II: Implementation and Tests in DUNE. Computing, 82(2–3):121–138, 2008a.
 Bastian et al. [2008b] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, M. Ohlberger, and O. Sander. A Generic Grid Interface for Parallel and Adaptive Scientific Computing. Part I: Abstract Framework. Computing, 82(2–3):103–119, 2008b.
 Bastian et al. [2010] P. Bastian, G. Buse, and O. Sander. Infrastructure for the coupling of Dune grids. In Proc. of ENUMATH 2009, pages 107–114. Springer, 2010.
 Bear [1988] J. Bear. Dynamics of Fluids in Porous Media. Dover Publications, 1988.
 Berninger et al. [2011] H. Berninger, R. Kornhuber, and O. Sander. Fast and robust numerical solution of the Richards equation in homogeneous soil. SINUM, 49(6):2576–2597, 2011.
 Borouchaki et al. [2000] H. Borouchaki, P. Laug, and P. George. Parametric surface meshing using a combined advancingfront – generalizedDelaunay approach. Int. J. Numer. Meth. Eng., 49(1–2):233–259, 2000.
 Bressan et al. [2014] A. Bressan, S. Čanić, M. Garavello, M. Herty, and B. Piccoli. Flows on networks: recent results and perspectives. EMS Surv. Math. Sci., 1(1):47–111, 2014.
 Cattaneo and Zunino [2014a] L. Cattaneo and P. Zunino. Computational models for fluid exchange between microcirculation and tissue interstitium. Networks and Heterogeneous Media, 9(1):135–159, 2014a.
 Cattaneo and Zunino [2014b] L. Cattaneo and P. Zunino. A computational model of drug delivery through microcirculation to compare different tumor treatments. International Journal for Numerical Methods in Biomedical Engineering, 30(11):1347–1371, 2014b.
 D’Apice et al. [2010] C. D’Apice, S. Göttlich, M. Herty, and B. Piccoli. Modeling, simulation, and optimization of supply chains. SIAM, 2010.
 Doussan et al. [1998] C. Doussan, L. Pages, and G. Vercambre. Modelling of the Hydraulic Architecture of Root Systems: An Integrated Approach to Water Absorption—Model Description. Annals of Botany, 81(2):213–223, 1998. doi: 10.1006/anbo.1997.0540.
 Dunbabin et al. [2013] V. M. Dunbabin, J. A. Postma, A. Schnepf, L. Pagès, M. Javaux, L. Wu, D. Leitner, Y. L. Chen, Z. Rengel, and A. J. Diggle. Modelling root–soil interactions using threedimensional models of root growth, architecture and function. Plant and Soil, 372:93–124, 2013. doi: 10.1007/s111040131769y.
 Dziuk and Elliott [2007a] G. Dziuk and C. M. Elliott. Finite elements on evolving surfaces. IMA J. Numer. Anal., 27:262–292, 2007a.
 Dziuk and Elliott [2007b] G. Dziuk and C. M. Elliott. Surface finite elements for parabolic equations. J. Comput. Math., 25(4):385–407, 2007b.
 Dziuk and Elliott [2013] G. Dziuk and C. M. Elliott. Finite element methods for surface PDEs. Acta Numerica, 22:289–396, 2013.
 D’Angelo [2007] C. D’Angelo. Multiscale modelling of metabolism and transport phenomena in living tissues. Bibliotheque de l’EPFL, Lausanne, 2007.
 Engwer and Müthing [to appear] C. Engwer and S. Müthing. Concepts for flexible parallel multidomain simulations. In Domain Decomposition Methods in Science and Engineering XXII. Springer, to appear.
 Flemisch et al. [2011] B. Flemisch, M. Darcis, K. Erbertseder, B. Faigle, A. Lauser, K. Mosthaf, P. Nuske, A. Tatomir, M. Wolff, and R. Helmig. DuMuX : DUNE for Multi{ Phase, Component, Scale, Physics, …} Flow and Transport in Porous Media. Advances in Water Resources, 34:1102–1112, 2011.
 Garavello and Piccoli [2006] M. Garavello and B. Piccoli. Traffic flow on networks. American Institute of Mathematical Sciences (AIMS), 2006.
 Geuzaine and Remacle [2015] C. Geuzaine and F. Remacle. Gmsh Reference Manual, 2015. www.geuz.org/gmsh/doc/texinfo/gmsh.pdf.
 Gräser and Sander [2009] C. Gräser and O. Sander. The dunesubgrid module and some applications. Computing, 8(4):269–290, 2009.
 Gross and Reusken [2011] S. Gross and A. Reusken. Numerical Methods for Twophase Incompressible Flows. Springer, 2011.
 Javaux et al. [2008] M. Javaux, T. Schröder, J. Vanderborght, and H. Vereecken. Use of a ThreeDimensional Detailed Modeling Approach for Predicting Root Water Uptake. Vadose Zone Journal, 7(3):1079, 2008. doi: 10.2136/vzj2007.0115.
 Kedem and Katchalsky [1958] O. Kedem and A. Katchalsky. Thermodynamic analysis of the permeability of biological membranes to nonelectrolytes. Biochimica et biophysica Acta, 27:229–246, 1958.
 Lang et al. [2011] S. Lang, V. J. Dercksen, B. Sakmann, and M. Oberlaender. Simulation of signal flow in 3d reconstructions of an anatomically realistic neural network in rat vibrissal cortex. Neural Networks, 24(9):998–1011, 2011.
 Leitner et al. [2010] D. Leitner, S. Klepsch, G. Bodner, and A. Schnepf. A dynamic root system growth model based on LSystems. Plant and Soil, 332:177–192, Jan. 2010. ISSN 0032079X. doi: 10.1007/s1110401002847.
 McClure et al. [2015] M. McClure, M. Babazadeh, S. Shiozawa, and J. Huang. Fully coupled hydromechanical simulation of hydraulic fracturing in threedimensional discrete fracture networks. In SPE Hydraulic Fracturing Technology Conference, 3–5 February, The Woodlands, Texas, USA, 2015. doi: http://dx.doi.org/10.2118/173354MS.
 McClure and Horne [2013] M. W. McClure and R. N. Horne. Discrete Fracture Network Modeling of Hydraulic Stimulation – Coupling Flow and Geomechanics. Springerbriefs in Earth Sciences. Springer, 2013.
 Motti et al. [1986] E. D. Motti, H.G. Imhof, and M. G. Yasargil. The terminal vascular bed in the superficial cortex of the rat: an SEM study of corrosion casts. Journal of Neurosurgery, 65(6):834–846, 1986.
 Nammi et al. [2010] S. Nammi, P. Myler, and G. Edwards. Finite element analysis of closedcell aluminium foam under quasistatic loading. Materials and Design, 31:712–722, 2010.
 Nitschke et al. [2012] I. Nitschke, A. Voigt, and J. Wensch. A finite element approach to incompressible twophase flow on manifolds. J. Fluid Mech., 708:418–438, 2012.
 Olshanskii et al. [2009] M. A. Olshanskii, A. Reusken, and J. Grande. A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal., 47(5):3339–3358, 2009.
 Pagès et al. [2004] L. Pagès, G. Vercambre, and J. Drouet. Root Typ: a generic model to depict and analyse the root system architecture. Plant and Soil, 258:103–119, 2004.
 Quarteroni and Formaggia [2003] A. Quarteroni and L. Formaggia. Mathematical Modelling and Numerical Simulation of the Cardiovascular System. In Modelling of Living Systems, Handbook of Numerical Analysis Series. EPFL, 2003.
 Reuther and Voigt [2015] S. Reuther and A. Voigt. The interplay of curvature and vortices in flow on curved surfaces. SIAM Multiscale Model. Simul., 13(2):632–643, 2015.
 Rognes et al. [2013] M. E. Rognes, D. A. Ham, C. J. Cotter, and A. T. T. McRae. Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2. Geosci. Model Dev., 6:2099–2119, 2013.
 Schmidt and Siebert [2005] A. Schmidt and K. G. Siebert. Design of Adaptive Finite Element Software. The Finite Element Toolbox ALBERTA, volume 42 of LNCSE. Springer, 2005. albertafem.de.
 Secomb et al. [2000] T. Secomb, R. Hsu, N. Beamer, and B. Coull. Theoretical simulation of oxygen transport to brain by networks of microvessels: effects of oxygen supply and demand on tissue hypoxia. Microcirculation, 7(4):237–247, 2000.
 Somma et al. [1998] F. Somma, J. W. Hopmans, and V. Clausnitzer. Transient threedimensional modeling of soil water and solute transport with simultaneous root growth, root water and nutrient uptake. Plant and Soil, 202:281–293, 1998.
 T. Witkowski [2012] A. V. T. Witkowski, R. Backofen. The influence of membrane bound proteins on phase separation and coarsening in cell membranes. Phys. Chem. Chem. Phys., 14:14509–14515, 2012.
 Tyree [1997] M. T. Tyree. The CohesionTension theory of sap ascent : current controversies. Journal of Experimental Botany, 48(315):1753–1765, 1997.
 Tyree and Zimmermann [2002] M. T. Tyree and M. H. Zimmermann. Xylem Structure and the Ascent of Sap, volume 12. Springer, 2002. ISBN 9783540433545. doi: 10.1016/03781127(85)900817.
 Vey and Voigt [2007] S. Vey and A. Voigt. AMDiS: adaptive multidimensional simulations. Comput. Visual. Sci., 10:57–67, 2007.
 Wolff [2013] M. Wolff. Multiscale modeling of twophase flow in porous media including capillary pressure effects. PhD thesis, Universität Stuttgart, 2013.