An AFC-stabilized implicit finite element method for partial differential equations on evolving-in-time surfaces

Accepted Manuscript An AFC-stabilized implicit finite element method for partial differential equations on evolving-in-time surfaces Andriy Sokolov, Ramzan Ali, Stefan Turek PII: DOI: Reference:

S0377-0427(15)00147-8 http://dx.doi.org/10.1016/j.cam.2015.03.002 CAM 10056

To appear in:

Journal of Computational and Applied Mathematics

Received date: 24 September 2014 Revised date: 28 February 2015 Please cite this article as: A. Sokolov, R. Ali, S. Turek, An AFC-stabilized implicit finite element method for partial differential equations on evolving-in-time surfaces, Journal of Computational and Applied Mathematics (2015), http://dx.doi.org/10.1016/j.cam.2015.03.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An AFC-stabilized implicit finite element method for partial differential equations on evolving-in-time surfaces Andriy Sokolov, Ramzan Ali, Stefan Turek Institut f¨ur Angewandte Mathematik, TU Dortmund Vogelpothsweg 87 44227 Dortmund, Germany

Abstract In this article we present a new implicit numerical scheme for reactiondiffusion-advection equations on an evolving in time hypersurface Γ(t). The partial differential equations are solved on a stationary quadrilateral, resp., hexahedral mesh. The zero level set of the time dependent indicator function φ(t) implicitly describes the position of Γ(t). The dominating convectivelike terms, which are due to the presence of chemotaxis, transport of the cell density and surface evolution may lead to the non-positiveness of a given numerical scheme and in such a way cause appearance of negative values and give rise of nonphysical oscillations in the numerical solution. The proposed finite element method is constructed to avoid this problem: implicit treatment of corresponding discrete terms in combination with the algebraic flux correction (AFC) techniques make it possible to obtain a sufficiently accurate solution for reaction-diffusion-advection PDEs on evolving surfaces. Keywords: level set, evolving surfaces, FEM, FCT, TVD, membrane, pattern-formation 1. Introduction Current applications in biomathematics demand fast, efficient and accurate solvers for systems of reaction-diffusion-advection PDEs. For example, Email addresses: [email protected] (Andriy Sokolov), [email protected] (Ramzan Ali), [email protected] (Stefan Turek)

Preprint submitted to SI: ACOMEN 2014

February 27, 2015

pattern formations and colouring of the animals coat [30, 31, 32, 51], dynamics of the yeast cell polarity [17, 39] and other processes are modeled by reaction-diffusion equations for which the instability of the Turing type causes formation of complex patterns. Chemotaxis models are widely used to describe bacteria/cell aggregation and formation processes [1, 29, 48, 49, 50], modeling of tumor invasion and metastasis processes before a proliferation dominated stage [3, 6, 7, 8, 9], modeling of vasculogenesis [2, 16, 40], etc. Here, besides diffusion and reaction one deals with the advection-like term which is due to specie-specie or specie-agent interactions. During the last years there are tendencies to extend these mathematical models to more complex, but on the other hand more physically reasonable systems by coupling them with PDEs, which are defined on an evolving-intime manifold Γ(t), see, e.g., [4, 5, 19] for reaction-diffusion patterns on growing surfaces, [22, 28] for chemotaxis processes on stationary and deforming manifolds, [17, 39] for protein dynamics on a cell’s membrane, etc. Mathematically, aforementioned physical models can be described by the following system of reaction-advection-diffusion equations ∂ci = Dic ∆ci − ∇ · [χci wic (c, ρ)ci ] + fi (c, ρ), in Ω, (1) ∂t ∂ ∗ ρj = Djρ ∆Γ(t) ρj − ∇Γ(t) · (χρj wjρ (c, ρ)ρj ) + gj (c, ρ), on Γ(t),(2) ∂t ∂∗ρ

where ∂t j is a time-derivative, which takes into account the evolution of Γ(t) and will be explained in the next section. Vector fields wc (·) and wρ (·) in equations (1) and (2) can describe the species-chemoattractant or speciesspecies interaction (e.g., w = ∇cj ) or some external velocity field. χc and χρ are (chemo-, resp., specie-) sensitivities, which can be nonlinear. f (·) and g(·) are some kinetic terms. Here, ci (x, t), i = 1 . . . n, are defined in the whole domain Ω and are solutions of (1). Unknown functions ρj (x, t), j = 1 . . . m live on the surface Γ(t) ⊂ Ω and are solutions of (2). The corresponding initial and boundary (if any) conditions for ci and ρj have to be provided. We adopt the notation by writing vector fields in bold letters, i.e., c = (c1 , . . . , cn )T . The next sections are dedicated to present an efficient and accurate technique for the solution of equations (2), with the examples included. We assume that the solution ρ of (2) can be (naturally) prolonged from Γ(t) to an ǫ-band Ωǫ (t), see Figure 1. The domain of interest or also the calculational domain is Ω = Ωin ∪ Ωout ∪ Γ. The position of Γ(t) is prescribed implicitly by the zero level set of the time-dependent level set function φ(t), those values are either given analytically or found by solving 2

the transport equation ∂φ + v · ∇φ = 0, (3) ∂t where v is the velocity of the surface, which can be prescribed analytically, or v can also be determined through forces acting on the surface Γ, e.g. chemo-attractive, chemo-repulsive forces, bending surface etc. By the proper choice of numbers of equations n and m, parameters Dic and Djρ , chemotaxis/advection-related functions χci and vic , as well as kinetic terms fi (·) and gj (·), one can obtain a large variety of models to be used in the biomathematics. In a series of papers [43, 44, 45] the authors constructed a robust and efficient numerical scheme for chemotaxis probW (t) lems in 2 and 3 spatial dimensions, in the W in case when n = 2 and m = 1. There, the G(t) W out Flux-Corrected Transport (FCT) or Total Variation Diminishing (TVD) stabilization techniques, Newton-like solvers and coupled as well as decoupled approaches Figure 1: Geometric illustration were analyzed, see also [42]. It was shown that the solver was able to deliver physically appropriate and accurate numerical solutions. Using the FCT method and the operator splitting technique one can extend the proposed framework to models of multiple interactive species in the presence of several chemoattractants/repellents, see, e.g., [20, 21]. In the paper [41] we constructed a numerical scheme for the equation of the type (2) coupled with a chemotaxis system and the surface Γ was considered to be stationary. In this article we will propose a numerical scheme for solution of the partial differential equation of type (2) on the time dependent surface with j = 1. This method treats implicitly the arising advective terms and boundary integrals. The resulting numerical scheme is equipped with a robust, positivitypreserving finite element based TVD/FCT-stabilization technique. For alternative stabilization of partial differential equations on evolving in time surfaces we would like to refer the interested to, e.g., works of Olshanskii et al. [11, 37]. The article is organized as follows: in section 2 we present the required theoretical derivations. Namely, in subsection 2.2 we explain the advective surface derivative and surface-related parameters, and in subsection 2.1 we describe the level set methodology and present the corresponding integraltheorems. Section 3 deals with the construction of an implicit numerical scheme, discusses the treatment of surface-advection and boundary-integral e

3

terms and discusses the corresponding stabilization techniques. In the following section 4 we discuss embedding of the FCT/TVD stabilization techniques into the implicit scheme. Then, in section 5 we present numerical results for reaction-diffusion partial differential equations in two spatial dimensions and examine the order of time- and space-convergence. Section 6 summarizes the characteristics of the proposed approach. 2. Theoretical results 2.1. Level set method We assume that Γ(t) ⊂ Ω ⊂ Rd+1 is a compact connected and oriented ddimensional hypersurface without boundary and that Γ is C 2 . Then, there exists a two times continuously differential level set function φ : T × Ω → R defined by   < 0 , if x is inside Γ(t), φ(t, x) = = 0 , if x ∈ Γ(t), (4)   > 0 , if x is outside Γ(t),

such that |∇φ| 6= 0. An outward normal to Γ(t) is n = (n1 , . . . , nd+1 )T = and the matrix PΓ =

(

δij −

d+1 X j=1

ni nj

)

ij

=I−

∇φ |∇φ|

(5)

∇φ ∇φ ⊗ |∇φ| |∇φ|

(6)

denotes the orthogonal projection onto the tangent space Tx Γ. Observe that if φ(·) is chosen as a signed distance function then |∇φ| = 1. For a scalar function ξ on Ω and a tangential vector field ξ on Γ one obtains ) ( d+1 X ∂ξ ∂ξ − ni nj , (7) ∇Γ ξ = PΓ ∇ξ = ∂xi j=1 ∂xj i ! d+1 d+1 X ∂ξi X ∂ξi − ∇Γ · ξ = ni nj . (8) ∂x ∂x i j j=1 i=1

We would like to mention that for the natural extension, i.e., when ξ is constant along ∇φ, one can omit Γ in spatial differentiation on the right-handsides of (7) and (8). The Laplace-Beltrami operator on Γ(t) with respect to the level set function φ can be written as ∆Γ ξ = ∇Γ · ∇Γ ξ = ∇Γ · PΓ ∇ξ. 4

(9)

The Eulerian mean curvature is defined through the level set function as H = −∇Γ · n = −∇Γ ·

∇φ . |∇φ|

In the level set formulation one perform calculation not only on Γ, but in some surrounding of Γ bulk Ωε of the width ε. For the sake of brevity and simplicity, we assume that all considered in this article values/fields can be naturally, i.e., as a constant in a ∇φ-direction, extended to any level set Γ(t)c = {x ∈ Ωε : φ(t, x) = c},

(10)

where c is some constant. In the following R we will omit c-notation and will simply write, e.g., n instead of nc or Ωε ∇Γ(t) P · ∇Γ(t) ϕ|∇φ| instead of R ∇Γc (t) P · ∇Γc (t) ϕ|∇φ|. Furthermore, for the sake of brevity and clarΩε ity, we will also assume that Ωε = Ω. While this assumption can be too restrictive for real 3D-applications, for our test purposes this assumption is justified. Interested readers are referred to works of Olshanskii et al. for the trace FEM [33] and/or for the space-time FEM [35, 36] which are high order accurate in space and time, stable, optimal with respect to d.o.f. and able to handle topological changes. 2.2. Surface derivatives and related derivations In this article we derive an implicit numerical scheme for the following parabolic chemotaxis-like equation on an evolving hypersurface Γ(t) ∂ ∗ρ = D∆Γ(t) ρ − ∇Γ(t) · (χρ∇Γ(t) c) + g(ρ) on Γ(t) × T, ∂t

(11)

which can be obtained from (2) by setting j = 1 and Dρ = D,

χρ = χ,

wρ = ∇Γ(t) c,

g(c, ρ) = g(ρ).

Here, the derivative

∂ ∗ρ = ∂t• ρ + ρ∇Γ(t) · v (12) ∂t is due to the evolution of Γ(t) and can be obtained by the Leibniz formula Z Z d ρ= ∂t• ρ + ρ∇Γ(t) · v. (13) dt Γ(t) Γ(t) 5

By ∂t• ρ = ∂t ρ+v ·∇ρ one denotes the advective surface material derivative. The surface velocity v = V n + vS can be decomposed into velocity components in the normal direction V n, with n to be a surface outward normal vector, and in the tangential direction vS . One can easily find that V (x, t) = −

φt (x, t) . |∇φ(x, t)|

Using the relations ∇Γ ·v = ∇Γ V ·n+V ∇Γ ·n+∇Γ ·vS = V ∇Γ ·n+∇Γ ·vS = −V H +∇Γ ·vS , ∂ρ + vS · ∇ρ, ∂n where H is a mean curvature, we can rewrite (11) as v · ∇ρ = V n · ∇ρ + vS · ∇ρ = V

∂ρ + ρ∇Γ(t) · vS = ∂n = D∆Γ(t) ρ − ∇Γ(t) · (χρ∇Γ(t) c) + g(ρ),

∂t ρ + vS · ∇ρ − V Hρ + V

(14)

or, in terms of the surface material derivative, as ∂t• ρ + ρ∇Γ · v = D∆Γ(t) ρ − ∇Γ(t) · (χρ∇Γ(t) c) + g(ρ).

(15)

3. Numerical scheme Our method is based on the method proposed by Dziuk and Elliott in [15]. This algorithm is the level set-based finite element scheme, where diffusive and R advective terms are treated implicitly but the surface-evolution term P v · ϕ|∇φ| treated explicitly. This method has some interesting propΩ erties: first, stationary finite elements save time of calculation and reduce the complexity of a solver. Secondly, the level set nature makes it possible to couple surface-defined PDEs with the domain-defined PDEs. Additional prescription of a transport equation for the level set function can allow to model a complex behavior of the evolving in time hypersurface Γ(t). And thirdly, the method does not require numerical evaluation of the curvature. Nevertheless, for more practical use this method has to be improved because of a few drawbacks. On the one hand, this method does not include any stabilization for advective-like terms, thus leaving opening for the numerical scheme to loose its positivity-preserving properties and leading to oscillatory solutions with negative values. On the other hand, this method does not include integral-boundary terms. Under some conditions the absence of corresponding discrete terms in the numerical scheme may lead to the loss of 6

accuracy of the solution near boundaries. This may also cause loss of positivity of the solution near boundaries follows by deterioration of the entire solution. To overcome these stumbling-blocks, we adopted a fully implicit FCT/TVD-stabilized finite element method inclusive boundary integrals. In the following we will briefly remind the derivation of the method by Dziuk and Elliott, especially those parts their our modifications are done, and then we show how it is possible to couple it with the FCT/TVD stabilization techniques to guarantee the positivity and smoothness of the solution. In the subsequent we will need the following lemmas, proofs can be found in [10, 15]. Lemma 1. (Coarea formula) ¯ → R be Lipschitz continuous and assume Let for each t ∈ [0, T ], φ(t, ·) : Ω that for each r ∈ (infΩ φ, supΩ φ) the level set Γr = {x|φ(x, ·) = r} is a ¯ → R is smooth d-dimensional hypersurface in Rd+1 . Suppose that η : Ω continuous and integrable. Then  Z supΩ φ Z Z η dr = η|∇φ|. (16) infΩ φ

Γr

Ω

Lemma 2. (Eulerian integration by parts) Assume that the following quantities exist. Then, for a scalar function η and a vector field Q we have Z Z Z ∇Γ(t) η|∇φ| = − ηHn|∇φ| + η(n∂Ω − n · n∂Ω n)|∇φ|, (17) Ω

Z

Ω

Ω

∇Γ(t) · Q|∇φ| = − Z −

Ω Z

Z

Ω

∂Ω

HQ · n|∇φ| +

∇Γ(t) · Qη|∇φ| + Ω

Q · nηH|∇φ| +

Z

ZΩ

Z

∂Ω

Q · (n∂Ω − n · n∂Ω n)|∇φ|,(18)

Q · ∇Γ(t) η|∇φ| =

∂Ω

Q · (n∂Ω − n · n∂Ω n)η|∇φ|,

(19)

where n∂Ω is an outward normal to ∂Ω. Lemma 3. (Implicit surface Leibniz formula) Let η be an arbitrary level set function defined on Ω such that the following quantities exist. Then Z Z Z d • η|∇φ| = (∂t η + η∇Γ(t) · v)|∇φ| − ηv · n∂Ω |∇φ|. (20) dt Ω Ω ∂Ω 7

Let us denote ϕ to be some test function, then using the coarea formula (16) we can write the variational formulation of (15): Z Z
• D∆Γ(t) ρ − ∇Γ(t) · (χρ∇Γ(t) c) ϕ|∇φ| (∂t ρ + ρ∇Γ · v) ϕ|∇φ| = Ω ZΩ + g(ρ)ϕ|∇φ|. (21) Ω

Due to the implicit surface Leibniz formula (20) the left hand side of (21) can be rewritten as Z Z d • (∂t ρ + ρ∇Γ(t) · v)ϕ|∇φ| = ρϕ|∇φ| dt Ω Ω Z Z − ρ∂t• ϕ|∇φ| + ρϕv · n∂Ω |∇φ|. (22) Ω

∂Ω

In general the boundary integral cannot be neglected. In Section (5) we show that its non-inclusion into the numerical scheme can cause kinks and oscillations in the solution near boundaries which grow and spread
throughout the domain as time evolves. Since D∇Γ(t) ρ − χρ∇Γ(t) c · n = 0, integration by parts in (19) gives Z (D∇Γ(t) ρ − χρ∇Γ(t) c) · ∇Γ(t) ϕ|∇φ| = ΩZ − ϕ∇Γ(t) · (D∇Γ(t) ρ − χρ∇Γ(t) c)|∇φ| ZΩ + (D∇Γ(t) ρ − χρ∇Γ(t) c) · n∂Ω ϕ|∇φ|. (23) ∂Ω

R The boundary integral ∂Ω (D∇Γ(t) ρ − χρ∇Γ(t) c) · n∂Ω ϕ|∇φ| on the right hand side of (23) cannot be neglected in general as well. Its numerical treatment though can be performed in a similar way to the boundary integral R ρϕv · n∂Ω |∇φ| from (22). For the brevity everywhere in this article we ∂Ω assume that the boundary ∂Ω is aligned with some level set Γr and therefore (D∇Γ(t) ρ − χρ∇Γ(t) c) · n∂Ω = 0. This assumption might be justified by the fact, that in real aplication a Γ-bulk, in which the surface PDE is considered/calculated, is much smaller than Ω and therefore almost never reaches its boundary ∂Ω. Applying (22) and (23) to (21), we obtain Z Z d ρϕ|∇φ| + (D∇Γ(t) ρ − χρ∇Γ(t) c) · ∇Γ(t) ϕ|∇φ| = dt Ω Z ZΩ Z • ρϕv · n∂Ω |∇φ| + g(ρ)ϕ|∇φ|. (24) ρ∂t ϕ|∇φ| − Ω

∂Ω

Ω

8

For the discretization in space we use time-independent (conforming) bilinear (in 2d) or trilinear (in 3d) finite elements with the corresponding space of test functions Qh = span{ϕ1 , . . . , ϕN }. Therefore (25)

ϕ˙ = v · ∇ϕ

and the semi-discretization problem for (24) reads: find P ∈ Qh such that Z Z d P ϕ|∇φ| + (D∇Γ(t) P − χρ∇Γ(t) c) · ∇Γ(t) ϕ|∇φ| = (26) dt ΩZ Z Ω Z P v · ∇ϕ|∇φ| − P ϕv · n∂Ω |∇φ| + g(ρ)ϕ|∇φ| ∀ϕ ∈ Qh . Ω

∂Ω

Ω

Numerical solution of (26) requires careful treatment of advective-like and boundary terms on its right- and left-hand sides. We propose the following implicit version of the algorithm. Fully-discrete scheme of the implicit version: Given P m at the tm -the time instance and the time step ∆t = tm+1 − tm , then solve for P m+1 Z Z 1 m+1 m+1 P ϕ|∇φ | + (D∇Γm+1 P m+1 − wP m+1 ) · ∇Γm+1 ϕ|∇φm+1 | ∆t Ω Ω Z Z m+1 m+1 m+1 − P v · ∇ϕ|∇φ |+ P m+1 ϕv m+1 · n∂Ω |∇φm+1 | Ω

=

1 ∆t

∂Ω

Z

Ω

P m ϕ|∇φm | +

Z

Ω

(27)

gϕ|∇φm |,

for all ϕ ∈ Sh . The discrete space Sh is generated by the nodal basis functions ϕi , i = 1, . . . , N , Sh = span{ϕ1 , . . . , ϕN }. w is some vector field which advects the solution, e.g., set w = χ∇Γ(t) c to obtain (26). Choosing {ϕi } test functions as quadrilateral, resp., hexahedral finite elements in two- or three-dimensional space, the matrix form of the equation (27) looks like follows: [M (|∇φm+1 |) + ∆tL(D|∇φm+1 |) − ∆tK(wm |∇φm+1 |) − ∆tN (v m |∇φm+1 |) + ∆tR(|∇φm+1 |)] P m+1 = M (|∇φm |)P m + ∆tg(|∇φm |) 9

(28)

Here, M (·) denotes the consistent mass matrix, L(·) is the discrete LaplaceBeltrami operator, K(·) is the discrete on-surface advection operator with the linearized velocity wm , N (·) is the discrete operator due to the surface evolution and R(·) are discrete boundary integrals with the entries defined by the formulae Z mij (ψ) = ϕi ϕj ψ, (29) ZΩ lij (ψ) = PΓ ∇ϕi · ∇ϕj ψ, (30) Ω Z ϕi ψ · PΓ ϕj , (31) kij (ψ) = ZΩ nij (ψ) = ϕi v m+1 · ∇ϕj ψ, (32) Ω Z rij (ψ) = ϕi ϕj v m+1 · n∂Ω ψ , (33) ∂Ω Z gi (ψ) = ϕi gi ψ. (34) Ω

It is known that under some conditions the discrete Laplace-Beltrami operator L(D|∇φ|) can deteriorate and leads to an ill-conditioned matrix. In this article we do not investigate this situation and refer to the corresponding work of Olshanskii and Reusken [34]. Also relevant for this work are the papers which modify the extended Laplace-Beltrami operator in a bulk domain to avoid degeneracy, see [12, 18, 38]. The considered by us implicit method (27) satisfies mass-conservation properties, for detailed derivations see [15]. In this context it is worth to mention other methods for partial differential equations of evolving surfaces which conserve the mass of the solution, see [12] and [14]. The convective-like matrices K(·) and N (·) depending on velocities w and v, as well as the boundary-integral operator R(·) may lead to the loss the positivity-preserving properties and thus lead to nonphysical oscillations in the solution profile. This fact necessitates the use of some stabilization based techniques to treat the operators or terms which caused oscillation. Rewriting the equation (28) as [M (|∇φm+1 |) + ∆tL(D|∇φm+1 |) − ∆tK(|∇φm+1 |)] P m+1 = M (|∇φm |)P m + ∆tg(|∇φm |)

(35)

we apply the linearized flux-corrected transport algorithm to the operator K(|∇φm+1 |) = K(wm |∇φm+1 |) + N (v m+1 |∇φm+1 |) − R(|∇φm+1 |). 10

(36)

4. Stabilization Here, we present a short description of the embedding of the FCT/TVD stabilized algorithm into our implicit algorithm (27), or more precisely into its discrete matrix-written form (35). As shown by Kuzmin et al. [24, 25, 26], positivity constraints can be readily enforced at the discrete level using a conservative manipulation of the matrices M = {mij } and K = {kij }, supposing that the source term ∆tg n does not course any threat to positivity. The former is approximated by its diagonal counterpart ML constructed using row-sum mass lumping X ML := diag{mi }, mi = mij (|∇φ|). (37) j

Next, all negative off-diagonal entries of K are eliminated by adding an artificial diffusion operator D. For conservation reasons, this matrix must be symmetric with zero row and column sums. For any pair of neighboring nodes i and j, the entry dij is defined as [24, 25] ( max{−kij , 0, −kji }, j 6= i, dij = (38) P − k6=i dik , j = i. It is clear that dij = dji . The result is a positivity-preserving discretization of low order. By construction, the difference f between the residual of this scheme and that of the underlying Galerkin approximation F = (ML − M (|∇φ|))

P n+1 − P n − D P n+1 ∆t

(39)

admits a conservative decomposition into a sum of skew-symmetric antidiffusive fluxes   X d Fi = fij , fij = mij + dij (Pi − Pj ) = −fji , ∀j 6= i. (40) dt j6=i To achieve high resolution while keeping the scheme positivity-preserving we will apply the limiting process. Its idea is that each flux is multiplied by a solution-dependent correction factor αij ∈ [0, 1] and is inserted into the right-hand side of the non-oscillatory low-order scheme. The original Galerkin discretization corresponds to the setting αij := 1. It may be used in regions where the numerical solution is smooth and well-resolved. The setting αij := 0 is appropriate in the neighborhood of steep fronts. One says that the flux limiting is of the TVD type, if the antidiffusive flux limiting takes into account only modification of the convective operator K, i.e. 11

F = −D P n+1 . One says that the flux limiting is of the FCT type, if the antidiffusive flux limiting takes into account not only modification of the convective operator K, but also the corresponding modification of the mass P n+1 − P n − D P n+1 . matrix M , i.e. F = (ML − M (|∇φ|)) ∆t In essence, the off-diagonal entries of the sparse matrices M and K are replaced by m∗ij := αij mij ,

kij∗ := kij + (1 − αij )dij ,

while the diagonal coefficients of the flux-corrected Galerkin operators are given by X X m∗ii := mi − αij mij , kii∗ := kii − (1 − αij )dij . j6=i

j6=i

In implicit FEM-FCT schemes [24, 25, 26], the optimal values of αij are determined using Zalesak’s algorithm [52]. The limiting process begins with canceling all fluxes that are diffusive in nature and tend to flatten the solution profiles. The required modification is: fij := 0

if

fij (Pj − Pi ) > 0,

where P is a positivity-preserving solution of low order [24, 25, 26]. The remaining fluxes are truly antidiffusive, and the computation of αij involves the following algorithmic steps: (a) Compute the sums of positive/negative antidiffusive fluxes into node i X X Si+ = max{0, fij }, Si− = min{0, fij }. j6=i

j6=i

(b) Compute the distance to a local extremum of the auxiliary solution P Q+ i = max{0, max(Pj − Pi )}, j6=i

Q− i = min{0, min(Pj − Pi )}. j6=i

(c) Compute the nodal correction factors for the net increment to node i     m i Q− m i Q+ i i − + , Ri = min 1, . Ri = min 1, ∆tSi+ ∆tSi− (d) Check the sign of the antidiffusive flux and apply the correction factor ( min{Ri+ , Rj− }, if fij > 0, αij = min{Ri− , Rj+ }, otherwise. 12

To summarize this section we present the flow chart of required algebraic manipulations of the FCT/TVD stabilization technique for the implicit fullydiscrete scheme (35)– (36): 1. Obtain the low-order scheme [ML + ∆tL(D|∇φm+1 |) − ∆tL(|∇φm+1 |)]P m+1 = ML P m + ∆tg(|∇φm |)

(41)

by performing mass lumping M → ML and transforming the highorder operator K into its nonoscillatory low-order counterpart L by adding a discrete diffusion operator D. 2. Remove excessive diffusion where possible by adding the limited antidiffusion operator F (α) [ML + ∆tL(D|∇φm+1 |) − ∆tL(|∇φm+1 |)]P m+1 + F (α) = ML P m + ∆tg(|∇φm |)

(42)

where F (α) is constructed according to the limiting process described above and α’s are found from algorithmic steps (a)–(d). In our test examples we will distinguish the following cases: 1. the implicit scheme without stabilization (IS) corresponds to (28). 2. the explicit scheme without stabilization (ES) corresponds to [M (|∇φm+1 |) + ∆tL(D|∇φm+1 |) − ∆tK(wm |∇φm+1 |) + ∆tR(|∇φm+1 |)] P m+1 = M (|∇φm |)P m + ∆tN (v m |∇φm |)P m + ∆tg(|∇φm |)

(43)

3. the implicit scheme without boundary term (ISwB) corresponds to (28) [M (|∇φm+1 |) + ∆tL(D|∇φm+1 |) − ∆tK(wm |∇φm+1 |) − ∆tN (v m+1 |∇φm+1 |)] P m+1 = M (|∇φm |)P m + ∆tg(|∇φm |) 13

(44)

4. Implicit scheme with TVD stabilization (IS-TVD) corresponds to (42), where F (α) is of the TVD-type. 5. the implicit scheme with FCT stabilization (IS-FCT) corresponds to (42), where F (α) is of the FCT-type. Due to the level-set methodology and the algebraic nature of the flux-corrected schemes, TVD and some FCT methods can be implemented for the equation (35) without their major changes (one should operate with proper M (·), L(·) and K(·) matrices). Nevertheless, there are certain FCT schemes which necessitate the evaluation of ρxi during the antidiffusive flux-limiting (42), e.g. Gradient-based slope limiting in [27]. In this case the corresponding projection onto the tangential space Tx Γ is required. Numerical tests showed that the TVD scheme has worse accuracy than the FCT algorithm: P n+1 − P n term results in a smeared soneglecting the (ML − M (|∇φ|)) ∆t lution. Therefore, the FCT-based algorithms are in this case more preferable. Positivity constraints are enforced using a nonlinear blend of high- and loworder approximations through an algebraic manipulation of the matrices M and K. The limiting strategy is fully multidimensional and applicable to (multi-)linear finite element discretizations on unstructured meshes. The approach was successfully tested for domain-defined equations for chemotaxis in 2- and 3-spatial dimensions [43, 44] and for chemotaxis equations which are defined on stationary surfaces [41]. For detailed presentation of the FEM-FCT methodology, including theoretical analysis (stability, positivity, mass preservation and convergence) and technical implementation details (data structures, matrix assembly) the interested readers are referred to [23, 24, 25, 26, 27] and other related publications of Kuzmin et al. 5. Numerical results 5.1. Example 1 In the first example we would like to demonstrate that numerical stabilization of certain PDEs on surfaces is necessary. For these purposes we consider the transport equation ∂ t ρ + v · ∇Γ ρ = 0

(45)

on the unit sphere Γ = {x : |x| = 1}. Applying the level set methodology, we solve the transport equation (45) on level sets Γc ⊂ Ω. For the sake of 14

simplicity our calculational domain Ω is chosen to be a union of all level sets Γc with c ∈ [0.5, 1.5]. The following initial condition ( 10 if |x − (0, 0, 1)T | ≤ 0.3 , ρ(x, t) = (46) 0 else. and the advective velocity vector-field v = {x1 , 0, −x3 }T

are taken. The mesh is constructed by refining the coarsest level via connecting opposite midpoints several times. In the table below we give the number of cells and degrees of freedom at every level of refinement. For level cells 1 192 2 1 536 3 12 288 4 98 304 786 432 5

d.o.f. 250 1 746 13 090 101 442 798 850

Table 1: Mesh levels on the computational domain Ω.

numerical experiments we use a grid of trilinear finite elements at the 5th level of mesh refinement. In Figure 2(a) we show the initial condition for ρ. Its numerical solution for the pure Galerkin scheme at an exemplary timepoint T = 0.2 is demonstrated in Figure 2(b). One can clearly observe artificial oscillations and negative values of ρ near regions of steep gradients. These nonphysical negative values grow rapidly as time evolves, which leads to the divergence of the algorithm and the corresponding termination of the simulation run. In [41] we also showed that the pure Galerkin scheme for chemotaxis problems on stationary surfaces cannot guarantee positivity preservation and smoothness of the solution. The corresponding FCT/TVD methodology helps to stabilize this type of problems and delivers a sufficiently accurate solution, see Figures 2(c) and 2(d)Ras well as [41]. In the following we show that the integral boundary term ∂Ω ρϕv · n∂Ω |∇φ| can also lead to undesired kinks near the boundary of the calculation domain, which can also spread through the whole domain and spoil the numerical solution. An advective term to the surface evolution requires the special numerical treatment as well. In the following example 5.2 we will show, that its pure Galerkin discretization without any stabilization method might lead to artificial oscillation in numerical solution, especially in regions with steep gradients. For corresponding reconstruction processes we refer, e.g., to [47]. 15

(a) initial solution

(b) pure Galerkin method

(c) TVD

(d) FCT

Figure 2: Numerical results for the transport problem, ∆t = 0.001.

5.2. Example 2 We test the following PDE on surfaces around a narrow band for pure Galerkin, FCT and TVD stabilization schemes ∂ ∗ρ + αρ = D∆Γ(t) ρ ∂t

on

Γ(t).

(47)

Using equation (12) and the Leibniz formula (13), one can obtain the following equation ∂ρ + v · ∇ρ + ρ∇Γ · v + αρ = D∆Γ(t) ρ ∂t

on Γ(t),

where α = 0.2. Let us suppose the initial prescription of the level set function as φ(x, t) = |x| − (1.0 + b t sin(5γ))

16

with b = 10 and γ = atan2(x2 , x1 ), defined as:  arctan( xx12 ) ,    arctan( x2 ) + π , x1 atan2(x2 , x1 ) = x2  arctan( )−π ,  x1   π , +2

x1 x2 x2 x2

> 0 ≥ 0, x1 < 0 < 0, x1 < 0 > 0, x1 = 0

(48)

(a) level set t = 0.00

(b) level set, t = 0.02

(c) level set, t = 0.05

(d) initial solution

(e) IS, lev. 5

(f) IS-TVD, lev. 5

(g) IS-FCT, lev. 5

(h) IS, lev. 6

(i) IS-TVD, lev. 6

(j) IS-FCT, lev. 6

(k) IS, lev. 7

(l) IS-FCT, lev. 7

Figure 3: Comparison of IS, IS-TVD and IS-FCT, red color stands for ρ = 1 and the blue for ρ = −1.

17

The computational domain is Ω = {x ∈ R2 : 0.5 ≤ |x| ≤ 1.5} and the initial condition ρ0 is ( 10.0 if x − (1.0.0) ≡ 0, ρ0 = 0.0 otherwise. The figures 3(a)-3(c) are plotted to show evolution of the level set at different time instances. Figure 3(d) is the initial solution and corresponding results in figures 3(e)-3(g) are obtained through standard Galerkin (IS scheme), ISTVD and IS-FCT stabilization schemes, respectively. Here, we use mesh at 5th level of refinement (which corresponds to 16640 d.o.f) and with a small diffusion coefficient D = 10−5 . Numerical oscillations in the solution profile of the Galerkin scheme are clearly seen here. At the same time TVD and FCT schemes deliver smooth and oscillatory-free solutions. Later, in figures 3(h)-3(j), we consider the same tests-setting but calculated on a finer refinement level 6 and with a larger diffusion coefficient D = 10−2 . It appears that the mesh refinement is still coarse for the standard Galerkin scheme to over come nonphysical peaks in the numerical profile. Finally, figures 3(k) and 3(l) represent the same case but with high diffusion coefficient D = 1.0. Also, in this case the mesh refinement is two times finer. One can observe that all schemes (IS, IS-TVD/FCT) deliver similarly smooth numerical solution, here we omit TVD result to avoid redundancy because it looks like the result obtained via FCT. In figures 3(e)-3(l), corresponding numerical results are calculated at time instance T = 0.05 with fixed time step ∆t = 0.0005. ¿From these test cases one concludes the convective term, which is due to evolution of the surface, might lead to oscillation in the vicinity of Γ-band, especially for steep gradients of ρ and coarse mesh refinements. Therefore some numerical stabilization is needed. The proposed AFC-based schemes of TVD and FCT types yield a smooth, positivity preserving and physically suitable numerical solution. 5.3. Example 3 Let us take example 1 from [15]. We solve the following equation ∂ ∗ ρ(x, t) = D∆Γ(t) ρ(x, t) + g(x, t) ∂t

on

Γ(t),

(49)

where Γ(t) is prescribed by the zero level set of the function φ(x, t) = |x| − 1.0 + sin(4 t)(|x| − 0.5)(1.5 − |x|).

(50)

As a domain we choose Ω = {x ∈ R2 : 0.5 ≤ |x| ≤ 1.5}. The boundary of the domain ∂Ω is aligned with a curve from the family Γr . The analytical 18

solution is chosen to be 2

ρ(x, t) = e−t/|x|

x1 . |x|

(51)

Since Γ(t) is time-dependent, the equation (49) transforms into ∂t ρ + vS · ∇ρ + V

∂ρ − V Hρ + ρ∇Γ · vS − ∆Γ ρ = g, ∂n

(52)

where H is a mean curvature of Γ(t) and therefore H = −1/|x|. Substituting vS = 0 into (52) we get ∂t ρ + V

∂ρ − V Hρ − ∆Γ ρ = g. ∂n

(53)

The function ρ(x, t) from (51) solves ∂t ρ − ∆Γ ρ = 0. Therefore, one finds that ∂ρ − V Hρ = V ρ g=V ∂n



1 2t + 3 |x| |x|



.

As the initial condition we set ρinit = ρ(x, t = 0). In figures 4(c) and 4(d) we demonstrate numerical results after the 100th iteration with the fixed time step ∆t = 0.0005 obtained by a scheme, which does not include a boundary integral term. Here, explicit and implicit schemes reveal the same artifact in the numerical solution: artificial kinks near the outside boundary of the domain might propagate in time to the whole domain causing deterioration of the numerical solution. In figures 5(a), 5(b) and 5(c) we show numerical solutions, which are obtained by the implicit scheme (27). One can clearly see much better profiles of the numerical solution. The 3rd level of the mesh refinement is used. In figure 5(a) we demonstrate the numerical solution of (27) without any stabilization. In figure 5(b) the TVD and in figure 5(c) the FCT stabilization techniques are applied. In table 2 we show the absolute errors and the obtained orders of convergence for the implicit scheme without stabilization 27, resp. 28. Since for this test case no numerical oscillations are produced, the embedded FCT and TVD schemes deliver almost same solutions as those of the pure Galerkin scheme. The corresponding errors are defined by (cf. [15]) L2 (Ωǫ )-error =



1 |Ωε |

Z

2

Ωǫ

(ρanalyt (x, t) − ρnum (x, t)) |∇φ| 19

 21

,

(a) mesh at level 2 and initial solution

(b) analytical solution no name

(c) numerical solution, level 3

(d) boundary kinks, level 2

Figure 4: The ISwB scheme causes kinks near boundaries. The 2nd and 3rd levels of mesh refinement consists of 288 and 1088 degrees of freedom respectively.

H1 (Ωǫ )-error =



1 |Ωε |

Z

2

Ωǫ

(∇Γ ρanalyt (x, t) − ∇Γ ρnum (x, t)) |∇φ|

20

 21

.

(a) IS scheme

(b) IS-TVD scheme

no name

(c) IS-FCT scheme

(d) no boundary kinks

Figure 5: Numerical solutions obtained by adoption of the fully implicit scheme (27).

lev.

d.o.f

2 3 4 5

288 1088 4224 16640

2 3 4 5

288 1088 4224 16640

2 3 4 5

288 1088 4224 16640

L2 (Ωǫ ) -errors order ε = 0.5 2.60532E-003 – 6.76148E-004 1.946 1.79701E-004 1.911 4.93317E-005 1.865 ε = 0.25 1.65152E-003 – 4.12974E-004 1.999 1.02201E-004 2.014 2.44469E-005 2.063 ε = 0.125 1.28775E-003 – 3.21906E-004 2.000 7.90935E-005 2.025 1.90064E-005 2.057

H1 (Ωǫ ) -errors

order

3.61167E-002 1.80981E-002 9.15163E-003 4.68961E-003

– 0.996 0.983 0.964

2.44530E-002 1.22096E-002 6.10269E-003 3.05108E-003

– 1.00 1.00 1.00

1.95126E-002 9.74873E-003 4.82267E-003 2.42365E-003

– 1.001 1.015 0.992

Table 2: L2 (Ωǫ ) and H1 (Ωǫ ) errors, order of convergence and number of degrees of freedom in Ωε .

21

5.4. Example 4 As a 4th test case we take example 2 from [15]: we solve the equation (49) in the domain Ω = {x ∈ R2 : 0.5 ≤ |x| ≤ 1.5} on stationary level sets φ(x, t) = |x| − 1.0. Here, the initial solution is ρ0 (x) = sin(4γ), where γ ∈ [0, 2π] is an angle between the x1 -axis and x = (x1 , x2 )T , defined in equation (48) and the tangential velocity of the surface Γ is vS = 0. Since γt = 0, the normal component of the surface velocity V is also zero. The mean value of ρ0 vanishes on every level set Γr , hence the solution tends to zero as time tends to infinity. But this occurs at a rate which depends on the radius of the circle because of the different diffusion coefficients on the different circles. Obtained numerical solutions at successive time instances are presented in figure 6. Here, IS, IS-FCT and IS-TVD schemes deliver same numerical results.

(a) at t = 0.00

(b) at t = 0.002

(c) at t = 0.045

(d) at t = 0.1

Figure 6: Solution of example 4 at various time instances, ∆t = 0.0001

22

5.5. Example 5 Now everything is similar to the previous example 5.4, but the tangential velocity of the surface is defined as vS = 10

(−φx2 , φx1 ) . |∇φ|

(54)

Numerical results at some instances of time intervals are shown in figure 7.

(a) at t = 0.0

(b) at t = 0.002

(c) at t = 0.045

(d) at t = 0.1

Figure 7: Solution of example 5 at various time instances, ∆t = 0.0001

In table 2 we show the L2 (Ωǫ ) and H1 (Ωǫ ) errors with order of convergence and number of degrees of freedom with different mesh refinement levels on the stripes 0.5 < |x| < 1.5, 0.75 < |x| < 1.25 and 0.85 < |x| < 1.15, respectively. 6. Conclusion Modern biological and medical applications require numerical solutions of systems of partial differential equations (PDEs) on manifolds/surfaces 23

and their coupling with domain-defined PDEs. The level set methodology largely allows for such a coupling by a surface-surrounding bulk and provides the mathematical machinery for the development of numerical methods for PDEs on closed manifolds. Substantial difficulties occur when the surface and its corresponding level set function are time dependent, i.e., non-stationary. Here, one does not only have to carefully treat additional terms which arise due to the surface evolution, but also has to come up with a robust stabilization for advective terms which are caused by deformations of the surface, chemotaxis effects, advective surface material derivatives, boundary integrals, etc. In this article we proposed a fully implicit numerical scheme for systems of reaction-diffusion-convection/chemotaxis equations on an evolving in time surface. The positivity preservation is guaranteed by the TVD/FCT stabilization technique. The presented numerical results show a sufficient accuracy, given a domain Ω ⊂ R2 . The employed mathematical techniques such as the level set method, finite element discretization in space and TVD/FCT stabilization of advective terms provide a straightforward extension of a given method to the 3-dimensional space with an embedded 2-dimensional hypersurface. Next steps will be to implement and comprehensively study the given method in 3D in order to verify its merits for more complex bio-medical applications. The proposed numerical scheme was implemented in the open source FEM-library FEAT2, which is developed and maintained at the department of Mathematics at the TU Dortmund. The downloading of this software and the corresponding tutoring is possible either from www.featflow.de or by writing an email to one of the authors of the article. Acknowledgments Ramzan Ali is supported through faculty development program from ”University of Central Asia” in collaboration with ”DAAD”. [1] M. Aida, T. Tsujikawa, M. Efendiev, A. Yagi and M. Mimura, Lower estimate of the attractor dimension for a chemotaxis growth system. J. London Math. Soc., 74, no. 2 (2006), 453–474. [2] D. Ambrosi, F. Bussolino and L. Preziosi, A review of vasculogenesis models. Computational and Mathematical Methods in Medicine: An Interdisciplinary Journal of Mathematical, Theoretical and Clinical Aspects of Medicine, 6, no. 1 (2005), 1–19. [3] A.R.A. Anderson, M.A.J. Chaplain, E.L. Newman, R.J.C. Steele and A.M. Thompson, Mathematical modelling of tumour invasion and metastasis. J. Theoretical Medicine 2 (1999), 129–154. 24

[4] R. Barreira, C.M. Elliott and A. Madzvamuse, The surface finite element method for pattern formation on evolving biological surfaces. J. Math. Bio., 63, no. 6 (2011), 1095–1119. [5] M. Bergdorf, I.F. Sbalzarini and P. Koumoutsakos, A Lagrangian particle method for reaction-diffusion systems on deforming surfaces. J. Math. Biol., 61, no. 5 (2009), 649–663. [6] M.A.J. Chaplain, The mathematical modelling of tumour angiogenesis and invasion. ACTA Biotheoretica, 43, no. 4 (1995), 387–402. [7] M.A.J. Chaplain, Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Bio., 60 (1998), 857–900. [8] M.A.J. Chaplain, Mathematical modelling of angiogenesis. Journal of Neuro-Oncology, 50 (2000), 37–51. [9] M.A.J. Chaplain and A.M. Stuart, A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor. IMA J. Math. App. Medicine and Bio., 10 (1993), 149–168. [10] I. Chavel, Riemannian Geometry: A Modern Introduction , Cambridge Studies in Advanced Mathematics, Cambridge University Press (1995). [11] A.Y. Chernyshenko and M.A. Olshanskii, An adaptive octree finite element method for PDEs posed on surfaces. arXiv:1408.3891, (2014). [12]

K. Deckelnick, C. Elliott and T. Ranner, Unfitted finite element methods using bulk meshes for surface partial differential equations, SIAM Journal on Numerical Analysis 52, no. 4 (2013), doi: 10.1137/130948641.

[13] B. Dong, Applications of Variational Models and Partial Differential Equations in Medical Image and Surface Processing. PhD theis, UCLA (2009). [14] G. Dziuk, C.M. Elliott, Finite elements on evolving surfaces. IMA Journal of Numerical Analysis, 27, no. 2 (2007), pp. 262–292, doi: 10.1093/imanum/drl023. [15] G. Dziuk, C.M. Elliott, An Eulerian approach to transport and diffusion on evolving implicit surfaces. Comput Visual Sci, 13, (2010), 17–28. 25

[16] A. Gamba, D. Ambrosi, A. Coniglio, A. de Candia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi, and F. Bussolino, Percolation, morphogenesis, and Burgers dynamics in blood vessels formation. Phys. Rev. Lett., 90 (2003). [17] A.B. Goryachev and A.V. Pokhilko, Dynamics of Cdc42 network embodies a Turing-type mechanism of yeast cell polarity. FEBS Lett 582, no. 10 (2008), 1437–1443. [18] J.B. Greer, An improvement of a recent Eulerian method for solving PDEs on general goemetries. J. Sci. Comput. 29, (2006), 216–246. [19] G. Hetzer, A. Madzvamuse and W. Shen, Characterization of Turing diffusion-driven instability on evolving domains. Discrete and Continuous Dynamical Systems - Series A, 32, no. 11 (2012), 3975–4000. [20] D. Horstmann, Generalizing the Keller-Segel Model: Lyapunov Functionals, Steady State Analysis, and Blow-Up Results for Multi-species Chemotaxis Models in the Presence of Attraction and Repulsion Between Competitive Interacting Species. J. Nonlinear Sci. 21 (2011), 231–270. [21]

D. Horstmann and M. Lucia, Nonlocal elliptic boundary value problems related to chemotactic movement of mobile species. RIMS, Kokyuroku Bessatsu B15 (2009), 39–72.

[22] C.M. Elliott, B. Stinner and C. Venkataraman, Modelling cell motility and chemotaxis with evolving surface finite elements. Journal of the Royal Society Interface, 9, no. 76,(2012), 3027–3044. [23] D. Kuzmin, R. L¨ohner and S. Turek, Flux-Corrected Transport, Springer, 2nd edition, 2012. [24] D. Kuzmin and M. M¨oller, Algebraic flux correction I. Scalar conservation laws, in “Flux-Corrected Transport: Principles, Algorithms, and Applications” (Eds. D. Kuzmin, R. L¨ohner, S. Turek), Springer, Berlin (2005), 155–206. [25] D. Kuzmin and S. Turek, Flux correction tools for finite elements. J. Comput. Phys., 175 (2002), 525–558. [26] D. Kuzmin, Explicit and implicit FEM-TVD algorithms with flux linearization. J. Comput. Phys., 228 (2009), 2517–2534.

26

[27] D. Kuzmin, Linearity-preserving flux correction and convergence acceleration for constrained Galerkin schemes. J. Comput. Appl. Math., 236 (2012), 2317–2337. [28] C. Landsberg, F. Stenger, A. Deutsch, M. Gelinsky, A. R¨osen-Wolff and A. Voigt, Chemotaxis of mesenchymal stem cells within 3D biomimetic scaffolds–a modeling approach, J. Biomech., 44, no. 2 (2011), 359–364. [29] M. Mimura and T. Tsujikawa, Aggregating pattern dynamics in a chemotaxis model including growth. Physica A, 230 (1996), 499–543. [30] J.D. Murray, Discussion: Turing’s theory of morphogenesisIts influence on modelling biological pattern and form. Bull. Math. Biol., 52 (1990), 119–152. [31] J.D. Murray, Mathematical Biology I: An Introduction. Springer Verlag, Third Edition (2002). [32] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications. Springer Verlag, (2003). [33] M.A. Olshanskii, A. Reusken and J. Grande, A Finite Element method for elliptic equations on surfaces. SIAM J. Numer. Anal. 47, no. 5 (2009), 3339-3358. [34] M.A. Olshanskii and A. Reusken, A finite element method for surface PDEs: Matrix properties. Numerische Mathematik 114, (2009), 491– 520. [35] M.A. Olshanskii, A. Reusken and X. Xu, An Eulerian space-time Finite Element method for diffusion problems on evolving surfaces. SIAM J. Numer. Anal. 52, no. 3 (2014), 1354-1377. [36] M.A. Olshanskii and A. Reusken, Error analysis of a space-time finite element method for solving PDEs on evolving surfaces. SIAM J. Numer. Anal. 52, no. 4 (2014), 2092-2120. [37] M.A. Olshanskii, A. Reusken and X. Xu, A stabilized finite element method for advection-diffusion equations on surfaces. IMA J. Numer. Anal. 28, (2014), 732–758. [38] M.A. Olshanskii and D. Safin, A narrow-band unfitted finite element method for elliptic PDEs posed on surfaces. submitted, (2014). 27

[39] A. R¨otz, M. R¨ager, Turing Instabilities in a Mathematical Model for Signaling Networks. J. Math. Bio., 65, 6-7 (2012), 1215–1244. [40] G. Serini, D. Ambrosi, E. Giraudo, A. Gamba, L. Preziosi and F. Bussolino, Modeling the early stages of vascular network assembly. The EMBO Journal, 22 (2003), 1771–1779. [41] A. Sokolov, R. Strehl and S. Turek, Numerical simulation of chemotaxis models on stationary surfaces. Discrete and Continuous Dynamical Systems, Series B, 18, no. 10 (2013), 2689–2704. [42] R. Strehl, Advanced numerical treatment of chemotaxis driven PDEs in mathematical biology. Phd thesis, TU Dortmund (2013). [43] R. Strehl, A. Sokolov, D. Kuzmin and S. Turek, A flux-corrected finite element method for chemotaxis problems. Computational methods in applied mathematics, 10, no. 2 (2010), 219–232. [44] R. Strehl, A. Sokolov, D. Kuzmin, D. Horstmann and S. Turek, A positivity-preserving finite element method for chemotaxis problems in 3D. Journal of Computational and Applied Mathematics, 239 (2013), 290–303. [45] R. Strehl, A. Sokolov and S. Turek, Efficient, accurate and flexible Finite Element solvers for Chemotaxis problems. Computers and Mathematics with Applications, 34, no. 3 (2011), 175–189. [46] L. Tian, C.B. Macdonald, and S.J. Ruuth, Segmentation on surfaces with the closest point method. , Int Conf on Image Processing, Cairo, Egypt, ICIP09 (2009). [47] S. Turek, O. Mierka, S. Hysing and D. Kuzmin, Numerical Study of a High Order 3D FEM-Level Set Approach for Immiscible Flow Simulation Numerical Methods for Differential Equations, Optimization, and Technological Problems. Comput. Methods Appl. Sci., 27 (2013), 65–91. [48] R. Tyson, S.R. Lubkin and J.D. Murray, A minimal mechanism for bacterial pattern formation. Proc. R. Soc. Lond. B, 266 (1999), 299– 304. [49] R. Tyson, S.R. Lubkin and J.D. Murray, Model and analysis of chemotactic bacterial patterns in a liquid medium. J. Math. Biol., 38 (1999), 359–375. 28

[50] R. Tyson, L. G. Stern and R. J. LeVeque, Fractional step methods applied to a chemotaxis model, J. Math. Biol., 41 (2000), 455–475. [51]

V.K. Vanag and I.R. Epstein, Pattern formation mechanisms in reaction-diffusion systems. Int. J. Dev. Biol. 53, no. (5-6) (2009), 673– 681.

[52] S.T. Z ALESAK, Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys., 31, (1979), 335–362.

29