Marc Goerigk
email: goerigk@mathematik.uni-kl.de

Sacha Hesse
email: sascha.hesse@student.uni-halle.de

Matthias Müller-Hannemann

Marie Schmidt
email: m.schmidt@math.uni-goettingen.de

Anita Schöbel
email: schoebel@math.uni-goettingen.de

Sascha Heße

Recoverable Robust Timetable Information

Keywords: 1998 ACM Subject Classification F.2.2 Nonnumerical Algorithms and Problems, G.2.2 Graph Theory (Graph algorithms, Network problems) timetable information, recoverable robustness, delay scenarios Digital Object Identifier 10.4230/OASIcs.ATMOS.2013

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In timetable information, the following problem is typically considered: Given a timetable, an origin and destination, and an earliest departure time, find the "best" route leading from origin to destination; see [START_REF] Müller-Hannemann | Timetable information: Models and algorithms[END_REF] for a survey. An obvious criterion to evaluate the quality of a route is its duration (or travel time); however, many other criteria have been suggested, as, e.g., the number of changes or the ticket costs [START_REF] Schnee | Fully realistic multi-criteria timetable information systems[END_REF][START_REF] Disser | Multi-criteria shortest paths in timedependent train networks[END_REF][START_REF] Berger | Fully dynamic speed-up techniques for multi-criteria shortest paths searches in time-dependent networks[END_REF]. Also the reliability of a path has been considered as a means to account for delays [START_REF] Disser | Multi-criteria shortest paths in timedependent train networks[END_REF][START_REF] Müller-Hannemann | Efficient timetable information in the presence of delays[END_REF][START_REF] Schnee | Fully realistic multi-criteria timetable information systems[END_REF]. In [START_REF] Dibbelt | Intriguingly simple and fast transit routing[END_REF], decision trees for passengers' travels under uncertainty are constructed. In a recent work [START_REF] Goerigk | The price of robustness in timetable information[END_REF][START_REF] Goerigk | The price of strict and light robustness in timetable information[END_REF], approaches from the field of robust optimization were considered. Robust optimization is an approach to handle uncertainty in optimization problems that dates back to the 70s [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF].

During the late 90s, it received new attention through the work of Ben-Tal, Nemirovski and co-authors [START_REF] Ben-Tal | Robust convex optimization[END_REF][START_REF] Ben-Tal | Robust solutions of linear programming problems contaminated with uncertain data[END_REF], that sparked a manifold of concepts and algorithms; among them the Γ-approach of [START_REF] Bertsimas | The price of robustness[END_REF], adjustable robustness [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF], light robustness [START_REF] Fischetti | Light robustness[END_REF], or recoverable robustness [START_REF] Liebchen | The concept of recoverable robustness, linear programming recovery, and railway applications[END_REF][START_REF] Stiller | Extending concepts of reliability. Network creation games, real-time scheduling, and robust optimization[END_REF]. In our work we focus on recoverable robustness. This is a two-stage concept: Given a set of recovery algorithms, a solution is considered as being robust when for every scenario it can be "repaired" using an recovery algorithm to become feasible. An application to the uncertain shortest path problem has been considered in [START_REF] Büsing | Recoverable robust shortest path problems[END_REF], where the set of recovery algorithms is given by exchanging up to a constant K arcs of the path. Related work can be found in [START_REF] Seref | Incremental network optimization: Theory and algorithms[END_REF], where a given path is updated to a new solution by either using or removing k arcs. Further applications of recoverable robustness include shunting [START_REF] Cicerone | Robust Algorithms and Price of Robustness in Shunting Problems[END_REF], timetabling [START_REF] Cicerone | Recoverable robustness in shunting and timetabling[END_REF][START_REF] Stiller | Extending concepts of reliability. Network creation games, real-time scheduling, and robust optimization[END_REF][START_REF] Liebchen | The concept of recoverable robustness, linear programming recovery, and railway applications[END_REF], platforming [START_REF] Caprara | Recoverable-robust platforming by network buffering[END_REF][START_REF] Stiller | Extending concepts of reliability. Network creation games, real-time scheduling, and robust optimization[END_REF][START_REF] Liebchen | The concept of recoverable robustness, linear programming recovery, and railway applications[END_REF], the empty repositioning problem [START_REF] Erera | Robust optimization for empty repositioning problems[END_REF], railway rolling stock planning [START_REF] Cacchiani | Railway rolling stock planning: Robustness against large disruptions[END_REF] and the knapsack problem [START_REF] Büsing | Recoverable robust knapsacks: the discrete scenario case[END_REF]. In some previous work [START_REF] Goerigk | The price of robustness in timetable information[END_REF][START_REF] Goerigk | The price of strict and light robustness in timetable information[END_REF], robust passenger information has been considered. It was shown that finding a strictly robust travel route which hedges against any possible delay scenario is an NP-hard problem and for practical application much too conservative. As an alternative, a robustness concept based on light robustness has been proposed. However, it is assumed that a passengers stays on the planned route whatever happens. In contrast to this, we allow that a passengers changes his/her route even if he/she already started the journey.

Contributions. In timetable information, as in many other problems, the passenger does not know the scenarios from the beginning of his/her trip, but learns the current scenario en route. This aspect has been neglected in previous work. In this paper, we describe a recoverable robustness approach to the timetable information problem which takes into account that the actual scenario is learned at some time point en route, and that the travel route may be updated from this point on. For such a recovery, all possible alternative routes may be chosen. The goal is to include this recovery step in the planning phase, i.e. to find a travel route which may be recovered for every delay scenario from a given uncertainty set.

Furthermore, our approach can deal with complicated delay scenarios, as they occur in public transportation where source delays cause the dropping of transfers and changes in the durations of driving and waiting activities. We develop polynomial-time algorithms that can handle any finite set of scenarios and test them on delay scenarios that are generated by propagating delay in transportation systems.

Using large-scale data modeling the train network of Germany, we show the effectiveness of our approach.

Overview. The remainder of this work is structured as follows: We shortly recapture the nominal timetable information problem, and introduce our recoverable robust model in Section 2. We present a polynomial-time label-setting algorithm in Section 3, and demonstrate its applicability to German railway data provided by Deutsche Bahn AG in Section 4. We conclude the paper and discuss further research directions in Section 5.

2

Model and Notation

Timetable information

In the following we refer to train timetables for the sake of simplicity; however, all results can be transferred to any other type of public transport. The starting point for our considerations is a directed acyclic graph, the so-called event-activity network (EAN) N = (E, A) which is regarded over a finite time horizon. Nodes E represent events in the train schedule: They can either be arrival events E arr (modeling the arrival of a certain train at a certain station), or departure events E dep (modeling the departure of a certain train from a certain station). Events are connected by directed arcs, the activities, which can be either driving activities A drive (modeling the trip of a train from a departure event to an arrival event), waiting activities A wait (modeling the time a train spends between an arrival and a departure event for passengers to embark and disembark), or transfer activities A trans (modeling passenger movements from one arrival event to another departure event within the same station). Each event i ∈ E has a schedule time π i ∈ IN; furthermore, to compute how delays spread within this network (see Section 4.2), we may assume that for each activity (i, j) ∈ A a minimal duration l ij is known, and thus a buffer time b ij := π j -π i -l ij . We assume that the initial timetable π is feasible, i.e., π j -π i ≥ l ij , hence all buffer times are nonnegative. The timetable information problem consists of finding a path within the event-activity network from one station to another, given an earliest departure time s. More precisely, we introduce two virtual events, namely one origin event u and one destination event v, corresponding to a given origin station s u and a destination station s v . The origin event u is connected by origin activities A org with all departure events at station s u taking place not earlier than s, while all arrival events at station s v are connected with v by destination activities A dest . We need to find a path P from u to v in N such that the nominal travel time t nom (P ) := π last(P ) -s on P is minimal, where last(P ) denotes the last arrival event on P .

Delays

Paths with minimal travel time in the EAN may be vulnerable to delays, i.e., in case of delays, the originally planned path may take much longer than planned, or planned transfers may even become infeasible if the connecting train does not wait for a delayed feeder train.

The aim of this paper is to give robust timetable information, i.e., to find paths in the EAN which are less vulnerable to delays. The delays observed in a public transportation system originate from source delays d a which can occur on the driving and waiting activities a of the train. These delays are partially absorbed by buffer times on the activities, however, they propagate through the network to subsequent events along driving and waiting activities and -if a transfer is maintained -along the corresponding transfer activity. We assume that each transfer is assigned a waiting time which specifies how long the connecting train will wait for the feeder train. If the delay of the feeder train exceeds the waiting time, the connecting train will depart on time. See Section 4.2 for details on our delay propagation method. We denote by A transfer (d) the set of maintained transfer activities in scenario d and denote the delay network N (d) := (E, A(d)) with A(d) := A drive ∪ A wait ∪ A transfer (d). The updated timetable is denoted by π(d). In this paper, we make the (simplifying) assumption that at some point in time, the passenger learns about all delays and can adapt ('recover') his/her travel route accordingly. We partition the events of the networks in a set U ξ of events where no delay has occurred so far and the passenger has not learned about future delays and a set V ξ where he/she knows all delays. We require the following properties of an information scenario ξ = (N ξ , π ξ , U ξ , V ξ ) consisting of a delay network N ξ , a disposition timetable π ξ on this network, and a partition (U ξ , V ξ ) of the events E:

u ∈ U ξ , v ∈ V ξ , if π ξ j > π j , j is in V ξ , all i with (i, v) ∈ A dest are contained in V ξ , if i is in V ξ , all successors of i are in V ξ .

AT M O S ' 1 3

A way to define the partition (U ξ , V ξ ) between nodes U ξ where no delay information is available and nodes V ξ with full delay information is to set U ξ := {j : π j < t ξ }, V ξ := {j : π j ≥ t ξ }, where t ξ denotes a revealing time t ξ ≤ min j∈E:π ξ j -πj >0 π j for every scenario ξ. For our computational experiments, we obtain N (ξ) := N (d ξ ) and π ξ := π(d ξ ) by delay propagation, see Section 4.2. However, our methods work for any set of scenarios ξ = (N ξ , π ξ , U ξ , V ξ ) as described above; it is not necessary to know the source delays to apply them. We define the set of activities where scenario ξ is revealed as A ξ := {(i, j) ∈ A : i ∈ U ξ , j ∈ V ξ }. A set of information scenarios will be called an uncertainty set and denoted by U. In this paper, we consider only finite uncertainty sets.

Recoverable Robust Timetable Information

Intuitively, we will call a path P recoverable robust if, when an information scenario ξ occurs while a passenger is traveling on P , this passenger can take a recovery path P ξ , to his/her destination. To formally define recoverable robust paths, we make use of the following observation: Let U be an uncertainty set and let P be a path from u to v in N .

Lemma 1.

For every ξ ∈ U, P contains exactly one arc from A ξ .

We denote this arc by (i ξ (P ), j ξ (P )). We denote by Q ξ (j) the set of recovery paths, i.e., all paths from a node j to v in N ξ , and set Q ξ (P ) := Q ξ (j ξ (P )). Definition 2. A path P is called recoverable robust (with respect to uncertainty set U) if for any ξ ∈ U the set of recovery paths Q ξ (P ) is not empty.

We assume that the passenger travels on the chosen path P until he/she learns about the information scenario he/she is in, i.e., until node j ξ (P ). Since at this node, the full information of ξ, i.e., N ξ , π ξ , U ξ and V ξ is revealed to the passenger, he/she can take the best path for this scenario. Thus, we assume that he/she reroutes from his/her current position according to scenario ξ.

The goal of this paper is to find "good" recoverable robust paths. However, there are different ideas on how to measure the quality of a recoverable robust path. We can evaluate the nominal quality: which recoverable robust path has shortest travel time if no delays occur? the worst-case quality: which recoverable robust path has the earliest guaranteed arrival time? Hence, we consider the following bicriteria problem: Problem 1. Bicriteria recoverable robust paths Input: EAN N = (E, A) with timetable π, origin u and destination v, starting time s, and uncertainty set U. Task: Find a path P from u to v in N which is recoverable robust and minimizes 1. the nominal travel time t nom (P ) = π last(P ) -s where last(P ) is the last arrival node on P (the nominal objective function) 2. the worst-case travel time t wc (P ) = max ξ∈U min Q∈Q ξ (P ) π ξ last(Q) -s where last(Q) is the last arrival event on Q (the worst-case objective function).

Note that for simplicity, we call t wc (P ) the worst-case travel time of P , although the path P is only taken in the nominal case and an alternative path P ξ := argmin Q∈Q ξ (P ) π ξ last(Q) -s is taken in case of delay scenario ξ.

In other words, the bicriteria recoverable robust shortest path problem aims at finding paths which, on the one hand, are good in the nominal case, i. e., if no delays occur, and on the other hand hedge against the scenarios from the uncertainty set U by minimizing worst-case travel time on the corresponding recovery paths.

3

Algorithms for Recoverable Robust Paths

A Recovery-Label Setting Algorithm

In this section we show that in case of finite uncertainty sets solutions to the bicriteria recoverable robust path problem can be found as solutions to a bicriteria minimax bottleneck shortest path problem in the EAN with recovery labels L(a) := (L nom (a), L wc (a)) T at all arcs a ∈ ξ∈U A ξ . The minimax bottleneck shortest path problem is the problem of finding a path between two nodes in a network which minimizes max a∈P c(a) in a graph with edge labels c(a). In the bicriteria version of this problem, every arc a is assigned two different labels c 1 (a) and c 2 (a). We now state the preprocessing Algorithm 1 which calculates the recovery labels L(a) := (L nom (a), L wc (a)) T needed to apply solution methods for the bicriteria minimax bottleneck shortest path problem. In Algorithm 1, for every a = (i, j) ∈ A ξ , L ξ (a) denotes the minimal travel time on a path which uses node j in scenario ξ. If no such path exists in scenario ξ, L ξ (a) is set to ∞. The algorithm returns the labels L = (L nom , L wc ) T which are 0 for all a / ∈ ξ∈U A ξ . For a ∈ ξ∈U A ξ , L nom (a) denotes the minimum nominal travel time when using a path containing node j, (and is ∞, if no such path exists) while L wc (a) represents the worst-case travel time for scenarios revealed at node j.

After initialization of all required labels to the value 0 (lines 1-6), we compute the shortest path distance from every event to the destination in the nominal scenario (line 7). This can be done by a single invocation of a standard shortest path tree computation in the reversed digraph from the destination v. Then, in the for-loop of lines 8-15, we iterate over all delay scenarios. With respect to the revealing time of scenario ξ, we now determine the set V ξ . Using again a backward shortest path tree computation with respect to N ξ , we determine for every event j ∈ E the length of a shortest path towards the destination v. Using these values, we can set the nominal and worst-case labels for paths which go through arcs in A ξ (lines 11-13). For ease of notation, we use ∞ + k = ∞ for all values k. Note that the label L nom (a) is only set if the corresponding edge a can be used in some scenario ξ ∈ U. We finally obtain the worst-case labels by taking the maximum over all scenarios. Note that lines 16-19 could be easily integrated into the main loop, but in the way presented here, the main loop can be run in parallel.

Given the recovery labels, the worst-case minimal travel time t wc (P ) on a path P can be calculated as the maximum over the labels L wc on P , as stated in the following lemma. Lemma 3. Let P be a path from u to v in N . Then for the labels calculated in Algorithm 1 it holds that if max a∈P L wc (a) < ∞, P is recoverable robust, and t wc (P ) = max a∈P L wc (a).

Proof. Consider an arbitrary scenario ξ := (N ξ , π ξ , U ξ , V ξ ). The passenger travels on path P until node j ξ (P ). Then, he/she can take the path calculated in step 10 of the algorithm until node v with total length L ξ (i ξ (P ), j ξ (P )) and this path has minimal length in N ξ among all paths containing node j. We conclude that (1) P is recoverable robust, and (2) t wc (P ) = max a∈P L wc (a).

For any path P , the labels L nom constitute lower bounds on the nominal travel time on P . However, for an arbitrary path P , the nominal traveling time can exceed max a∈P L nom (a). This can be avoided for paths which do not make detours after the scenarios are revealed.

AT M O S ' 1 3

Algorithm 1 Construction of recovery labels

Require: EAN N = (E, A) with timetable π, origin node u, destination node v, starting time s, and finite uncertainty set U. Ensure: Label L(a) ∈ R 2 + for every a ∈ A. 1: for (i, j) ∈ A do Initialization 2:

Set L nom (i, j) := 0.

3:

for ξ ∈ U do 4:
Set L ξ (i, j) := 0.

5:

end for 6: end for 7: Find length K nom (j) of shortest path from every j ∈ E to v in N . Set K nom (j) := ∞ if no such path exists. 8: for ξ ∈ U do 9:

Determine A ξ .

10:

Find length K ξ wc (j) of shortest path from every j ∈ E to v in N ξ . Set K ξ wc (j) := ∞ if no such path exists.

11:

for (i, j) ∈ A ξ do 12:

Set L nom (i, j) := π j -s + K nom (j).

Setting nominal labels.

13:

Set L ξ (i, j) := π ξ j -s + K ξ wc (j). Setting worst-case labels.

14:

end for 15: end for 16: for (i, j) ∈ A do 17:

Set L wc (i, j) := max ξ∈U L ξ (i, j)

18:

Set L(i, j) := (L nom (i, j), L wc (i, j)) T . 19: end for 20: return L Lemma 4. Let P be a path from u to v in N such that the path P 2 defined as the subpath of path P starting in the last arc (i, j) in P ∩ ξ∈U A ξ is a shortest path from j to v. Then for the labels calculated in Algorithm 1 it holds that if max a∈P L wc (a) < ∞, P is recoverable robust, t nom (P ) = max a∈P L nom (a), and t wc (P ) = max a∈P L wc (a).

Proof. This follows from Lemma 3 and from the construction of the labels L nom in Algorithm 1 as the sum of the travel time π j -s until node j and the shortest path travel time K nom (j) from j to v.

As a conclusion, we obtain the following theorem.

Theorem 5. The bicriteria recoverable robust path problem corresponds to a bicriteria bottleneck shortest path problem in the EAN with labels L.

It is folklore that the single-criteria bottleneck shortest path problem can be solved in linear time on directed acyclic graphs. The Pareto front of bicriteria bottleneck shortest path problems can be found in O(|A| 2 ) by a simple ε-constraint method which enumerates all possible values of the first objective function, deletes edges whose labels exceed the given value, and finds a bottleneck shortest path with respect to the second criterion in the remaining graph (compare [START_REF] De Lima Pinto | The tricriterion shortest path problem with at least two bottleneck objective functions[END_REF]).

Lemma 6. Algorithm 1 determines the labels L in time O(|A| • |U |).

Proof. The initialization takes time O(|A|•|U|).

Since we can assume that N is topologically sorted, shortest paths from a node to all other nodes can be found in time O(|A|). Hence, step 7 takes time O(|A|). For every ξ ∈ U, determining A ξ is in O(|A|). Since step 10 again is a shortest path calculation in a topologically sorted network and the operations in the loop over all (i, j) ∈ ξ∈U A ξ take constant time, steps 8-15 can be executed in time

O(|A| • |U |).

Single-Criteria Versions of Recoverable Robustness

To calculate the Pareto front of the bicriteria recoverable robust path problem with finite uncertainty set, we can use the approach as sketched in the previous section. However, we are also interested in two single-criteria versions of the problem. In particular, results of versions with single objective values can be much easier compared for sets of instances.

Problem 2. Worst-case optimal recoverable robust paths

Find a recoverable robust path P from u to v in N such that the nominal quality of P is smaller or equal than a given nominal quality bound T nom , P minimizes t wc (P ).

Problem 3. Nominally optimal recoverable robust paths Find a recoverable robust path P from u to v in N such that the worst-case quality of P is smaller or equal than a given worst-case quality bound T wc P minimizes t nom (P ).

Algorithm 2 describes how to compute worst-case optimal recoverable robust paths. The pseudo-code for an analogous algorithm to compute nominally optimal recoverable robust path, Algorithm 3, is provided in the Appendix.

Algorithm 2 Worst-case optimal recoverable robust path

Require: Network N = (E, A), labels L, nominal quality bound T nom , origin event u, destination event v. Ensure: Path P which is optimal for Problem 2 (if existing).

1: for a ∈ A do 2:

if L nom (a) > T nom then

3:

Remove a from A.

4:

end if 5: end for 6: Find a bottleneck shortest path P wc in N according to labels L wc . 7: if there is no such path with length < ∞ then 8:

return There is no recoverable robust path. 9: else 10:

Let (i, j) be the last arc on P wc ∩ ξ∈U A ξ .

11:

Denote by P 1 (j) the path P wc until node j.

12:

Find a shortest path P 2 (j) in N from j to v.

13:

return P := P 1 (j) ∪ P 2 (j), t nom (P ) := max a∈P L nom (a), t wc (P ) := max a∈P L wc (a)

14: end if

AT M O S ' 1 3

4

Experimental Results

Test Instances

The basis for our computational study is the German train schedule of February 1, 2013 from which we created an event-activity network. We generated transfer activities between pairs of trains at the same station provided that the departing train is scheduled to depart not later than 60 minutes after the planned arrival time of the feeding train. In addition, since some train lines operate only every two hours or irregularly, we add further transfer arcs. Namely, for each arrival event at some station s, we also create a transfer arc to those departure events which exceed the time bound of 60 minutes but provide the very next opportunity to get to a neighboring station. The main characteristics of the resulting network are shown in Table 1. To study the robustness of passenger paths, queries should not be too easy. For example, we are not interested in paths which do not require any transfer. Therefore, we decided to generate 1000 relatively difficult queries as follows. For each query, origin and destination are chosen uniformly at random from a set of the 3549 most important stations in Germany (this choice of stations has been provided by Deutsche Bahn AG). Such a pair of origin and destination stations is only accepted if the air distance between them is at least 200km and if the shortest travel route between them requires at least one transfer. The desired start time is uniquely set to 8:00am. The resulting set of queries has an average nominal travel time of 398 minutes and 3.3 transfers per query.

Generating Information Scenarios

A delay scenario d ∈ IN A drive ∪Await 0 specifies a delay on each driving and waiting activity. To generate a delay scenario, we first choose the revealing time of the scenario. Afterwards, we decide for each driving and waiting activity whether it shall receive a source delay or not. We use a parameter p ∈ (0, 1) specifying the probability that a train receives a source delay. This parameter p can be chosen depending on the level of robustness one wants to achieve.

If a train shall be source-delayed, we select one of its driving or waiting activities uniformly at random from those which are scheduled after the revealing time of the scenario and choose the source delay for this activity uniformly at random among 10, 15, 20, 25, and 30 minutes. The source delays on all other activities are set to 0. For simplicity, we assume that trains receive source delays independently from each other.

We use the following basic delay propagation rule in order to compute how delays spread along driving, waiting and maintained transfer activities: π(d) denotes the timetable adapted to delay scenario d. If the start event of an activity a = (i, j) is delayed, also its end event j will be delayed, where the delay can be reduced by the slack time b a . I.e. we require π(d) ≥ π and

π j (d) ≥ π i (d) + l a + d a (1)
for all activities a = (i, j) ∈ A wait ∪ A drive . For transfer activities equation ( 1) does not necessarily hold. Motivated by real-world decision systems of rail operators, we assume that the decision whether a transfer is actively maintained or not is specified by a fixed waiting time rule: Given a number wt a ∈ IN for every transfer activity, the transfer is actively maintained if the departing train has to wait at most wt a minutes compared to its original schedule. If transfer a is actively maintained, we require that (1) holds for it. However, if for a transfer activity a = (i, j) (1) holds due to some earlier delay on the train corresponding to j, a is maintained, even if π j (d) -π j > wt a . Hence, every delay d induces a new set of transfer activities which is denoted as A transf er (d). Given these waiting time rules for a given delay scenario d we can propagate the delay through the network along the activities in A drive ∪ A wait ∪ A trans (d) and, thus, calculate the corresponding adapted timetable according to the following propagation rule:

π j (d) = max π j , max i:(i,j)∈A; πi(d)+lij ≤πj +wtij {π i (d) + l ij + d ij } (2)
where we set wt a = ∞ ∀a ∈ A wait ∪ A drive and d a = 0 ∀a ∈ A trans . The concrete waiting time rule used in our experiments is that high speed trains (like Intercity Express ICE, Intercity IC, and Eurocity EC) wait for each other at most three minutes, whereas trains of other train categories do not wait. Note that delay propagation can be done in time O(|A|).

The uncertainty sets used in our experiments contain a number k of independent scenarios generated as described above.

Environment

All experiments were run on a PC (Intel(R) Xeon(R), 2.93GHz, 4MB cache, 47GB main memory under Ubuntu Linux version 12.04 LTS). Only one core has been used by our program. Our code is written in C++ and has been compiled with g++ 4.6.3 and compile option -O3.

AT M O S ' 1 3

Experiments

The purpose of this study is to evaluate the potential of recoverable robust paths as an alternative timetable information method in pretrip planning. A standard way of doing timetable information is to search for a path with minimum travel time as primary objective and with minimum transfers as a secondary one. We take this kind of standard search as the baseline of our comparisons.

Experiment 1: What is the effect of delays on the paths of the standard search?

We perform the following evaluation. Suppose that P is a given path. For each delay scenario, we determine the first event after the scenario's revealing time. We assume that the passenger can adjust his/her path to the delay scenario at this point and therefore compute the earliest arrival time at the destination under these conditions. The worst-case arrival time over all scenarios is the value we are interested in. To each of our 1000 test queries we applied the same set of 100 delay scenarios with parameter p = 0.20. We observe that on average the worst-case travel time is 450 minutes, i.e., 13% larger than the planned one. The absolute difference is 52 minutes on average.

Experiment 2: What is the price of a worst-case optimal recoverable robust path in comparison with a standard path? Using the same 100 delay scenarios as for Experiment 1, we are interested in two quantities, namely the nominal travel time and the worst-case travel time of a worst-case recoverable robust path. We upper bounded the nominal arrival time of a recoverable robust path by 150% of the fastest nominal path. Among all paths satisfying this bound we minimized the worst-case arrival time over all scenarios. Our computational results show that for all 1000 queries but two cases there exists a recoverable robust path. An interesting observation is that 34.2% of all standard paths are already the worst-case optimal recoverable robust paths. However, in 27% of the queries the worst-case arrival time is improved in comparison with the standard path. If there is an improvement, the reduction is 29 minutes on average, but the maximum observed difference is 220 minutes. The histogram in Figure 1 gives a more detailed picture. It shows how often a saving of x minutes in the worst-case scenario can be achieved by choosing a recoverable robust path. The price a passenger has to pay if he/she chooses a recoverable robust path is a slight average increase in nominal travel time to 407 minutes, i.e., about just 9 minutes more than for the standard search.

In Figure 2, we show box-and-whisker plots for the distributions of travel times for five algorithmic variants. The data is based on our test set of 1000 queries, each evaluated for 100 delay scenarios generated with parameter p = 0.2 for the probability that a train will be delayed by a source delay.

Recall that StNom and StWC stand for the nominal and worst-case travel time in minutes of the standard search, while RRNom and RRWC denote the nominal and worst-case travel time for worst-case optimal recoverable robust paths, respectively. Finally, SRNom gives the nominal travel time for strictly robust paths.

Experiment 3: What is the influence of parameter p, initially chosen as p = 0.2? Recall that parameter p specifies the probability that a train will be delayed by a source delay. To quantify the sensitivity of the different solution methods on the chosen uncertainty set, we redo the previous two experiments with p = 0.1 and p = 0.15. Figure 3 (left) and Table 2 summarize our findings and show the average nominal (StNom) and worst case travel time (StWC) in minutes for the standard search and the nominal (RRNom) and worst-case time (RRWC) for the optimal recoverable robust paths, respectively. If the probability parameter p increases, we observe a slight increase of average worst case travel times (what This histogram shows the number of cases where with respect to the worstcase scenario we can save x minutes by choosing a worst-case optimal recoverable robust path instead of the standard path. q q q q q q q q q q q q q q q q q q q q q q q StNom StWC RRNom RRWC SRNom should be expected), whereas the nominal travel time of recoverable robust paths is almost unchanged. We conclude that p = 0.2 might be preferable since it provides recoverability for the more severe scenarios at no price with respect to nominal travel time.

Table 2 shows the raw data from which Figure 3 (left) has been derived. Experiment 4: Comparison with strictly robust paths. Using the same uncertainty set as in the previous experiments, we computed the set of transfer activities which break at least once. We marked these arcs as forbidden, and rerun shortest path queries on the resulting even-activity network. Paths in this network are considered as strictly robust since no transfer will ever break. The average nominal travel time if we look for the fastest strictly robust path (SRNom) is 451 minutes for the uncertainty set with p = 0.2 (see also Figure 3 (left) and the last row of Table 2. Hence, the average nominal travel time of these paths is not better than the average worst-case time for standard paths. In full agreement with previous studies [START_REF] Goerigk | The price of robustness in timetable information[END_REF][START_REF] Goerigk | The price of strict and light robustness in timetable information[END_REF], strictly robust paths turn out to be too conservative. The parameter k = |U| has been chosen as a pragmatic compromise between efficiency (the computational effort scales linearly with k) and the degree of robustness we want to guarantee. Obviously, the more different scenarios we use, the higher the level of robustness we can achieve. Therefore, we fixed the parameter p = 0.20 but varied k ∈ {75, 100, 125}. Table 3 shows the average travel times in minutes for these variants, and Figure 3 (right) displays the additional travel time over the baseline of the standard path in minutes. It is interesting to observe that the average worst-case travel times depend only marginally on the parameter k in the chosen range. As expected, there is a slight increase of a few minutes on worst-case travel time when we increase k. At the same time, the average nominal travel time for recoverable robust paths does not increase. Further experiments will be needed to see whether this trend will be confirmed if k is chosen in an even wider range.

AT M O S ' 1 3

Practicality of our approach. For the purpose of this study, we have merely implemented a first prototype without much emphasis on performance issues. Our running times are several minutes per query which is clearly impractical. The main bottleneck is the computation of labels which grows linearly with the number of used scenarios. However, the most expensive part, namely the loop of lines 8-15, could be run in parallel. Thus, using massive parallelization and further speed-up techniques, we see a clear perspective that the computation time for a recoverable robust path can be brought down to a few seconds.

Conclusion and Further Research

In this work we introduced the concept of time-dependent recoverable-robust paths within the framework of timetable information. We showed that the resulting bicriteria problem can be solved in polynomial time using a label-setting algorithm, and a subsequent bottleneck shortest path calculation. The proposed concept and algorithm was experimentally evaluated on timetable information instances covering the whole German train network (schedule of 2013). While computation times are still too high for practical applications in the current implementation, we may assume that a parallelized algorithm will be sufficiently fast; moreover, as our experiments show that the proposed model has a valuable trade-off between nominal and worst-case travel times, such an algorithm will provide a customer-friendly alternative in practice. Further research includes the comparison of recoverable robust paths to lightly robust paths (see [START_REF] Goerigk | The price of robustness in timetable information[END_REF][START_REF] Goerigk | The price of strict and light robustness in timetable information[END_REF]), and the extension of the proposed model to multi-stage robustness where only partial information on the scenario is given at discrete points in time. Also, the evaluation of the computed paths with respect to a set of real delay scenarios is currently being analyzed.

A Algorithmic Approach

Algorithm 3 describes how to find nominally optimal recoverable robust paths subject to an upper worst-case quality bound.

Algorithm 3 Nominally optimal recoverable robust path

Require: Network N = (E, A), labels L, worst-case quality bound T wc , origin event u, destination event v. Ensure: Path P which is optimal for Problem 3 (if existing).

1: for a ∈ A do 2:

if L wc (a) > T wc then 3:

Remove a from A.

4:

end if 5: end for 6: Find a bottleneck shortest path P nom in N according to labels L nom . 7: if there is no such path with length < ∞ then 8:

return There is no recoverable robust path. 9: else 10:

Let (i, j) be the last arc on P nom ∩ ξ∈U A ξ .

11:

Denote by P 1 (j) the path P nom until node j.

12:

Find a shortest path P 2 (j) in N from j to v.

13:

return P := P 1 (j) ∪ P 2 (j), t nom (P ) := max a∈P L nom (a), t wc (P ) := max a∈P L wc (a) 14: end if

a∈PL

  wc (a) ≥ max a∈P L wc (a) = t wc (P ) (or t nom (P ) = max a∈P L nom (a) ≥ max a∈P L nom (a) = t nom (P ), respectively).

Figure 1

 1 Figure 1This histogram shows the number of cases where with respect to the worstcase scenario we can save x minutes by choosing a worst-case optimal recoverable robust path instead of the standard path.

Figure 2 Figure 3

 23 Figure 2Box-and-whisker plots for the travel time distributions of several algorithmic variants.

Table 1

 1 Characteristics of the used event activity network and test queries.

	characteristic	event activity network
	# trains	38,495
	# events	2,015,664
	# stations	8,857
	# transfer activities	19,869,867
	aver. nominal travel time	398 min
	aver. # transfers per query	3.3
	Lemma 7. Algorithm 2 and Algorithm 3 are correct.
	Proof. Let P be the path returned by Algorithm 2 or Algorithm 3. Then, due to the
	construction of P in step 13 of each algorithm, the assumptions of Lemma 4 are fulfilled,
	i.e.,
	since max

a∈P L wc (a) < ∞, P is recoverable robust, t nom (P ) = max a∈P L nom (a), and t wc (P ) = max a∈P L wc (a). Since there is no arc a with L nom > T nom (or L wc > T wc , respectively) we have that t nom (P ) ≤ T nom (or ( t wc (P ) ≤ T wc , respectively). Furthermore, for any other path P we have that t wc (P ) = max

Table 2

 2 Comparison of standard and robust solutions: Average travel time in minutes for k = |U | = 100 scenarios.

		p = 0.10 p = 0.15 p = 0.20
	StNom	398	398	398
	StWC	441	447	450
	RRNom	407	406	407
	RRWC	433	438	442
	SRNom	440	446	451

Table 3

 3 Comparison of standard and robust solutions for different sizes k of the uncertainty set: Average travel time in minutes for p = 0.20.

		k = 75 k = 100 k = 125
	StNom	398	398	398
	StWC	447	450	451
	RRNom	410	407	407
	RRWC	438	442	443

To which extent do our observations depend on the size of the scenario set?

  All previous experiments have been run with 100 different delay scenarios.

* partially supported by grants SCHO 1140/3-1 and MU 1482/4-3 within the DFG programme SPP 1307 Algorithm Engineering and by the European Union Seventh Framework Programme (FP7-PEOPLE-2009-IRSES) under grant number 246647 and the New Zealand Government (project OptALI). † Effort sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant number FA8655-13-1-3066. The U.S Government is authorized to reproduce and distribute reprints for Governmental purpose notwithstanding any copyright notation thereon.


P Díez

R Cottereau

S Zlotnik

A stable XFEM formulation for multi-phase problems enforcing the accuracy of the fluxes through Lagrange multipliers

When applied to diffusion problems in a multiphase setup, the popular XFEM strategy suffers from an inaccurate representation of the local fluxes in the vicinity of the interface. The XFEM enrichment improves the global quality of the solution but it is not enforcing any local feature to the fluxes. Thus, the resulting numerical fluxes in the vicinity of the interface are not realistic, in particular when conductivity ratios between the different phases are very high. This paper introduces an additional restriction to the XFEM formulation aiming at properly reproducing the features of the local fluxes in the transition zone. This restriction is implemented through Lagrange multipliers, and the stability of the resulting mixed formulation is tested satisfactorily through the Chapelle-Bathe numerical procedure. Several examples are presented and the solutions obtained show a spectacular improvement with respect to the standard XFEM.

INTRODUCTION

Multiphase problems appear in about all fields of physics and mechanics. The classical approach to solve these problems with the Finite Element method (FEM) consists in constructing a mesh that follows the interface. Therefore, each element of the mesh pertains to only one phase, and the solution of the coupled problem can be represented reasonably well by a polynomial approximation over each element. However, when the interface is moving, the constraint that the mesh should follow the interface means that the mesh must be reconstructed at each evolution step, which can rapidly become very costly. This problem is central in metal forging, oceanography, imaging, flame modeling, melting of materials, among other applications.

Similar concerns appear for single-phase problems with intricate geometries that are solved using the immersed boundary method [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF][START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Peskin | The immersed boundary method[END_REF], or similar techniques. In these approaches, a virtual interface is created where the boundary of the domain lied, while giving very stiff or soft properties to the newly-introduced material, depending on the type of boundary condition desired. This virtual two-phase problem is equivalent to the original, single-phase, problem, but the meshing constraints are relaxed. This set of methods is becoming very interesting, concurrently with the ever-widening use of real microscopy images for the definition of the geometry of computational problems [START_REF] Vennat | Collagen fiber network infiltration: permeability and capillary infiltration[END_REF].

In the two sets of problems described above (evolving multi-phase problems and geometrically intricate single-phase problems), the possibility that the elements of the mesh be intersected by the interface is very appealing. Using the FEM, it is technically possible to do so [START_REF] Barrett | Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces[END_REF][START_REF] Mackinnon | Treatment of material discontinuity in finite element computations[END_REF]. However, the rate of convergence of the solution with respect to the size of the elements is heavily deteriorated [START_REF] Pehlivanov | Superconvergence analysis of approximate boundary-flux calculations[END_REF]. Authors have therefore proposed alternative approaches. Among those, the Generalized and eXtended Finite Element Methods (GFEM/XFEM) have been widely developed [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF][START_REF] Oden | A new cloud-based hp finite element method[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Chessa | The extended finite element method (XFEM) for solidification problems[END_REF][START_REF] Merle | Solving thermal and phase change problems with the eXtended finite element method[END_REF][START_REF] Belytschko | Structured extended finite element methods of solids defined by implicit surfaces[END_REF][START_REF] Chessa | An extended finite element method for two-phase fluids: flow simulation and modeling[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF][START_REF] Dolbow | A numerical strategy for investigating the kinetic response of stimulus-responsive hydrogels[END_REF][START_REF] Duddu | A combined extended finite element and level set method for biofilm growth[END_REF][START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF][START_REF] Zlotnik | Hierarchical X-FEM for n-phase flow[END_REF][START_REF] Fries | The intrinsic XFEM for two-fluid flows[END_REF][START_REF] Cottereau | Modeling, with a unified level-set representation, of the expansion of a hollow in the ground under different physical phenomena[END_REF]. In these methods, the classical FEM basis, using polynomials over each element, is enriched with functions that incorporate information about the interface. In the particular problems considered here, the enrichment functions introduce the possibility for the solution to have a discontinuous gradient over the interface. The improvement of the XFEM solution is often dramatic in terms of global errors, with respect to a FEM solution. However, large errors in the evaluation of the fluxes close to the interface can arise, in particular when the conductivity ratio between the two phases is large. Unfortunately, these fluxes are often very important in practice. In particular, they often provide the main drive for the evolution of the interface [START_REF] Knapen | Resistance of soils to concentrated flow erosion: a review[END_REF]. The main objective of this paper is to propose a method to improve this flux evaluation in the vicinity of the interface, while retaining the advantages of using unfitted meshes.

Other authors have considered similar questions. In [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF], the authors considered enriched functions that consisted of the restriction of the classical finite element functions to each side of the interface. As these functions are not continuous over the interface, a variant of Nitsche's approach was used to weakly enforce the continuity of the primal variable. The main drawback of this formulation is that it involves a parameter that must be chosen with care in order to obtain a stable formulation. Nevertheless, several papers [START_REF] Harari | Analysis of an efficient finite element method for embedded interface problems[END_REF][START_REF] Zunino | An unfitted interface penalty method for the numerical approximation of contrast problems[END_REF] elaborated on variants of this original paper. Alternatively, some authors tried to enforce the continuity of the displacements through Lagrange multipliers [START_REF] Ji | On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method[END_REF]. The drawback of that approach is that it introduces additional variables in the elements cut by the interface. Also, care must be taken in the choice of the space of Lagrange multipliers in order to enforce stability. Several choices of Lagrange multiplier spaces have been proposed [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF], with the corresponding stability tested through numerical procedures [START_REF] Chapelle | The inf-sup test[END_REF][START_REF] Bathe | The inf-sup condition and its evaluation for mixed finite element methods[END_REF]. A mortar-like Lagrange approach was also proposed [START_REF] Kim | A mortared finite element method for frictional contact on arbitrary interfaces[END_REF] to allow for an independent discretization of the bulk and interface fields. Note that, to the best of our knowledge, all authors use enrichment functions that are discontinuous in both the displacement and the gradient. The enrichment functions that we use in this paper are discontinuous in the gradient but continuous in the displacement (see section 3).

Although it is not central to the issues addressed in this paper, we would lastly like to comment on the parameterization of the interface. Several methods have been proposed [START_REF] Crank | Free and moving boundary problems[END_REF], among which front tracking methods [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] and the marker-in-cell method [START_REF] Gorczyk | Large-scale rigid-body rotation in the mantle wedge and its implications for seismic tomography[END_REF][START_REF] Jan | A cell-by-cell thermally driven mushy cell tracking algorithm for phase-change problems[END_REF]. We consider here level-set functions [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Osher | Level set methods: an overview and some recent results[END_REF], which provide a very efficient and elegant alternative to the previous parameterizations. Since their first use in the description of dynamical two-phase fluid systems [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], their power has been acknowledged for the parameterization of complex evolving phases. In particular, their ability to deal with changes in topology without any remeshing has been recognized [START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF][START_REF] Karlsen | A fast marching method for reservoir simulation[END_REF]. They have been used in several fields of geophysics and geomechanics, including modeling of two-phase flows and permeability estimation in reservoir simulations [START_REF] Karlsen | A fast marching method for reservoir simulation[END_REF][START_REF] Nielsen | Reservoir description using a binary level set model[END_REF], tectonic plates subduction [START_REF] Zlotnik | Numerical modelling of tectonic plates subduction using X-FEM[END_REF][START_REF] Zlotnik | Modelling gravitational instabilities: slab break-off and Rayleigh-Taylor diapirism[END_REF], seismic waves travel time computation [START_REF] Sethian | 3-D traveltime computation using the fast marching method[END_REF], and, generally, for inverse problems and optimal design [START_REF] Ito | Level-set function approach to an inverse interface problem[END_REF][START_REF] Burger | A level set method for inverse problems[END_REF][START_REF] Burger | A survey of level set methods for inverse problems and optimal design[END_REF].

As stated earlier, the objective of this paper is to propose a method to improve the flux evaluation in the vicinity of the interface, while retaining the advantages of using unfitted meshes in the context of the XFEM. This paper builds on [START_REF] Cordero | XFEM+: una modificación de XFEM para mejorar la precisión de los flujos locales en problemas de difusión con conductividades muy distintas[END_REF], where the continuity of the gradients across the interface was enforced in a strong form. The resulting formulation was unfortunately very dependent on the type of elements that were used, and no stability analysis was possible. The general form of the mixed weak formulation presented at section 5 was then briefly sketched in [START_REF] Cottereau | Numerical modeling of erosion using an improvement of the extended finite element method[END_REF], but without the stability analysis and examples presented here. The adopted methodology is based in enforcing continuity of the flux across the interface using Lagrange multipliers. Consequently, the number of unknowns is increased with the dimension of the Lagrange multipliers space. Being the interface in a manifold of lesser order with respect to the computational domain, the cost increase is considered to be very moderate. In any case, the payoff in local accuracy worths the computational effort.

The outline of the paper is the following: in section 2, the problem of interest is stated; in section 3, the XFEM and the enrichment functions that we use are discussed; in section 4, an illustration of the lack of accuracy of the XFEM for the fluxes in the vicinity of the interface is presented; in section 5 and 6, which constitute the core of the paper, our Lagrange-based mixed formulation is introduced and its stability is discussed; and finally, in section 7, several examples illustrate the behavior of the method proposed. perature distribution over a bi-phasic material: find u such that:

Ω 1 Ω 2 Γ (a) P Γ n τ Ω 2 Ω 1 (b) Ω 1 Ω 2 Ω e 1 2 3 Γ int n (c)
-∇ • q = f in Ω 1 ∪ Ω 2 (1a) q = ν∇u in Ω 1 ∪ Ω 2 (1b) q • n = g N on Γ N (1c) u = u D on Γ D (1d) q • n = 0 on Γ (1e) u = 0 on Γ (1f)
In these equations, ν = ν 1 in the phase Ω 1 and ν = ν 2 in the phase Ω 2 . The normal vector n is the outgoing normal vector along Γ N ∪ Γ D and is oriented (arbitrarily) as indicated on figure 1(b), from Ω 1 towards Ω 2 . Note that the continuity of the (normal) flux, enforced by equation (1e), and the fact that the material coefficients ν 1 and ν 2 are different, implies that the gradient of u is necessarily discontinuous (in the direction of the normal to Γ), that is

∇u| Ω1 • n = ∇u| Ω2 • n. (2)
Note also that the jump of the normal component of the gradient depends on the contrast between ν 1 and ν 2 . On the contrary, the tangential component of the gradient is continuous, due to the continuity of u in Ω. Thus, denoting by τ the unit vector tangent to Γ at P (see figure 1(b)), one has

∇u| Ω1 • τ = ∇u| Ω2 • τ . (3)

Variational form

The following functional spaces are introduced to properly state the weak form of the problem. The space V containing the solution u is defined as

V := {u ∈ H 1 (Ω) : u = u D in Γ D }, (4)
and the corresponding test functions space, V 0 is

V 0 := {u ∈ H 1 (Ω) : u = 0 en Γ D }. (5)
Thus, the weak form of the problem, equivalent to (1) reads:

find u ∈ V such that a(u, v) = ℓ(v), for all v ∈ V 0 , (6)
where the bilinear and linear forms a(•, •) and ℓ(•) are given by

a(u, v) := Ω1∪Ω2 ν∇u • ∇v dΩ and ℓ(v) := Ω f v dΩ - Γ N g N v ds. ( 7
)
It is worth noting that the continuity constraint (1e) is implicitly imposed in [START_REF] Mackinnon | Treatment of material discontinuity in finite element computations[END_REF]. In fact, this condition is enforced in [START_REF] Mackinnon | Treatment of material discontinuity in finite element computations[END_REF] in a weak fashion, in the same manner as the Neumann boundary conditions (1c). Therefore, this restriction is exactly fulfilled only if the equation is solved exactly, using the infinitedimensional spaces V and V 0 . In the FEM solution, using finite-dimensional spaces approximating V H and V H0 , these conditions are only verified approximately. We come back to this issue in section 5.

PHASE TRACKING WITH LEVEL SETS AND XFEM ENRICHMENT

In this section, we briefly recall the definitions of level sets, of the FEM, and of the XFEM, with a particular emphasis on the type of enrichment functions that are used in our version of the XFEM.

Level sets

Level set functions [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Osher | Level set methods: an overview and some recent results[END_REF] provide a very efficient and elegant parameterization of multi-phases domains. In the simplest setting, they allow to discriminate between two areas of a domain, with no explicit parameterization of the actual interface. Conceptually, they are constructed in a space of higher dimension than the interface they intend to represent, with a smoother topological behavior that allows for an easier manipulation. The level set is a function φ(x) defined over the entire domain Ω, and its sign indicates the belonging to one or the other of the two areas Ω 1 and Ω 2 :

φ(x) =    > 0, x ∈ Ω 1 = 0, x ∈ Γ < 0, x ∈ Ω 2 (8)
Usually, its absolute value is defined as the distance to the interface. The interface is hence parameterized by the set of zeros of the level set.

Finite Element Method (FEM)

A finite element mesh generates a discrete functional space V H ⊂ V , where H stands for the characteristic size of the elements in the mesh. The mesh is a partition of the domain Ω into disjoint elements Ω e , e = 1, . . . , n e , that is Ω = e Ω e and Ω e ∩ Ω e ′ = ∅ for e = e ′ . The corresponding discrete counterpart of the test space V 0 is denoted by V H0 . The standard FEM basis of shape functions, generating V H0 is denoted by N i , i = 1, . . . , n H , being n H the number of nodal points in the finite element mesh. The affine space V H is therefore V H = {u ⋆ D } + V H0 , where u ⋆ D is a function in span{N 1 , . . . , N poin }, defined over Ω and that coincides with u D on Γ D . Typically, u ⋆ D is determined as the interpolation of u D with the FEM mesh. In the following, it is assumed that the FEM discretization properly reproduces u D (no oscillation terms in the data are considered). Thus, the FEM solution u H is defined by the nodal values u i , i = 1, . . . , n H :

u H = n H i=1 N i u i , (9)
and is such that

a(u H , v) = ℓ(v), for all v ∈ V H0 ⊂ V 0 . (10)
The discretization results in an algebraic system of equations

K uu u = F u , (11)
where the unknown column vector u contains the coefficients u i , associated to the nodes that are not on the Dirichlet boundary Γ D , and the stiffness matrix K uu and force vector

F u have generic components [K uu ] ij = a(N i , N j ) and [F u ] j = ℓ(N j ) -a(u ⋆ D , N j ).
A consequence of the assumption of no data oscillations is that the error associated with u H , e := uu H , belongs to V 0 . Thus, Galerkin orthogonality holds and u H is optimal in the sense that

u -u H ≤ u -v ∀ v ∈ V H ( 12
)
where • is the energy norm, defined by

v 2 = a(v, v).
The level set function φ is also approximated by φ H ∈ V H and determined by the nodal values φ

i , i = 1, . . . , n H , φ ≈ φ H = n H i=1 N i φ i . ( 13
)
This allows to describe the interface as a continuous line (surface in 3D), smooth inside the elements and with slope discontinuities when intersecting the element edges (sides in 3D). For linear triangular elements the interface is a polygonal line. Thus, the normal vector to the interface (see Figure 1(c)) is easily computed at any point of the interface in the interior of an element Ω e :

n = ∇φ H ∇φ H with ∇φ H = n H i=1 φ i ∇N i . (14)
Note that, as φ H is positive on Ω 1 , its gradient, and consequently n, point from Ω 1 towards Ω 2 .

EXtended Finite Element Method (XFEM)

The standard FEM provides a solution u H that is infinitely regular (smooth) inside any element Ω e . Therefore, standard finite elements are unable to reproduce the gradient jumps of the actual solution (see (1e) and ( 2)). The XFEM enriches the FEM solution using a ridge function, hence enabling the discretization to introduce a jump of the gradient across the interface. The ridge function used here is the one introduced in [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] and defined as

R = n H i=1 N i |φ i | - n H i=1 N i φ i . (15)
The function R vanishes in all the elements that are not crossed by the interface. Thus, the support of R is precisely the set of elements that should be enriched. In the following, the set of indices of these elements is denoted by E a . Introducing N a as the set of nodal indices corresponding to the nodes belonging to the enriched elements, the previous expression is rewritten as

R = i∈Na N i |φ i | - i∈Na N i φ i . (16)
For the sake of a simpler notation and without any loss in generality, it is assumed that the enriched nodes are the first in the node list, that is N a = {1, 2, . . . , card N a }. Thus, the XFEM approximation, u X , is

u X = n H i=1 N i u i + j∈Na RN j a j , (17)
where the coefficients a j , j ∈ N a , stand for the enriched nodal coefficients. The approximation u X lies in the discrete functional space V X = V H ⊕ span{RN i , i ∈ N a } and fulfills an equation analogous to [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] but replacing V H by V X . The resulting algebraic linear system of equations reads

K uu K ua K T ua K aa u a = F u F a , ( 18
)
where a is the column vector of unknowns a j , j ∈ N a and the generic entries of the matrices and vectors are [

K ua ] ij = a(N i , RN j ), [K aa ] ij = a(RN i , RN j ), and [F a ] j = ℓ(RN j ).
Similarly to the FEM solution u H , u X is optimal from an energetic point of view in the sense that Galerkin orthogonality stands in V X . Consequently, the error uu X is such that

u -u X ≤ u -v , ∀ v ∈ V X . ( 19
) Note that, taking v = u H ∈ V H ⊂ V X yields u -u X ≤ u -u H . ( 20
)

FLUX APPROXIMATION ON THE INTER-FACE

As previously stated, the FEM is unable to fulfill the flux continuity (1e) at the points on the interface in the interior of the elements. This is because the gradient of u H is continuous inside the elements and therefore equation ( 2) cannot be fulfilled at these points. Recall that the continuity of the normal flux at a point P of the interface reads

(q 2 • n 2 + q 1 • n 1 ) = (q 2 -q 1 ) • n = (ν 2 ∇u H -ν 1 ∇u H ) • n = 0. ( 21
)
However, the gradient of the standard finite element approximation, ∇u H , is continuous inside the elements, that is

(ν 2 ∇u H -ν 1 ∇u H ) • n = (ν 2 -ν 1 )∇u H • n. (22)
Thus, the continuity of the flux cannot be achieved with standard finite elements unless the gradient (and the flux) is identically zero. This is not the case in a general situation. Further, it should be noted that the approximation deviates especially far from the exact solution when the contrast between the two coefficients ν 1 and ν 2 increases. The XFEM uses an enriched basis that allows jumps in the normal component of the gradient of the solution ∇u H . Hence, it is expected that it may give u H the additional freedom to verify the flux continuity [START_REF] Duddu | A combined extended finite element and level set method for biofilm growth[END_REF]. Following equation [START_REF] Duddu | A combined extended finite element and level set method for biofilm growth[END_REF], the gradient of the approximation can be written:

∇u X = n H i=1 u i ∇N i + j∈Na (R∇N j + N j ∇R)a j . (23)
The possible discontinuity of ∇u X in a point P of Γ, in the interior of an element Ω e , appears thanks to the term ∇R. Indeed, ∇N i and R are both continuous inside Ω e . Further, since φ has been assumed positive inside the domain Ω 1 and negative inside Ω 2 (by convention, see equation ( 8)), and thanks to the definition of R in equation ( 16), we have the following expressions for R on both sides of the interface:

R| Ω1 = i∈Na N i (|φ i | -φ i ) and R| Ω2 = i∈Na N i (|φ i | + φ i ) . (24)
Note that ∇R| Ω2 -∇R| Ω1 = 2∇φ H = 0, where φ H denotes the level set function, as defined in equation [START_REF] Belytschko | Structured extended finite element methods of solids defined by implicit surfaces[END_REF]. As expected, the jump in the gradient of R is parallel to the normal n, since the latter is parallel to ∇φ H . Thanks to this additional freedom with respect to the classical finite element approximation, it is expected that the XFEM approximation may perform substantially better than the FEM.

In the rest of this section, we illustrate the behavior of the XFEM in terms of global error and quality of the approximation of the fluxes in the vicinity of the interface. We therefore present two simple examples: a square with a horizontal interface (intersecting the elements of the mesh), and the same with a tilted interface. In the first example, the XFEM captures the exact solution while the FEM behaves poorly. In the second example, the XFEM behaves well in terms of global error, but the fluxes are wrongly evaluated.

Exemple 1a (horizontal interface): Let us consider the problem (1), with the following parameters

: Ω = [-1, 1] × [-1, 1], g N = 0 for x = ±1, u D = 1 for y = 1, u D = 0 for y = -1, f = 0, ν 1 =
1000 and ν 2 = 1. Note that the ratio of conductivities is very large. We first consider a perfectly horizontal interface, at y = 0.1 (cutting the elements of the proposed mesh), with Ω 2 at y > 0.1. In this first configuration, the problem is quasi-1D. The exact solution is:

u(x, y) = 1 1100, 9 1000 (y + 1)
for y ≤ 0.1

y + 1099, 9 for y > 0.1 . ( 25
)
Since the solution is linear, if the mesh were to follow the interface, the FEM solution would be the exact one. However, in the case that we are interested in, the interface goes through the elements, and we can observe on figure 2 that the solution in terms of fluxes is rather inaccurate. More precisely, because the gradient of the FEM solution is continuous and the contrast in conductivities is very large, the fluxes appear to almost cancel on one side of the interface in the elements that are cut. Hence it appears that there are two jumps in the fluxes: one over the interface, created purely by jump in conductivities; and one between the line of elements cut by the interface and the next one, only on the weak side (where ν = 1). On the other side, this discontinuity between the elements cut by the interface and the first line perfectly included inside a phase is not so obvious. These last two remarks are related to the fact that, because it is not using a functional basis that incorporates information on the interface, the FEM effectively sees an average value of the conductivity rather than two phases. As the contrast if very large, this average conductivity is close to ν 1 = 1000, even for elements almost entirely contained in Ω 2 . Finally, it is important to understand that this local quantity of interest will not necessarily be better estimated with a more refined mesh. The global approximation will become better, but very locally, these two discontinuities will still be observed.

On this first configuration, on the other hand, the XFEM captures the analytical solution. Indeed, the enriched basis includes the exact solution. Note that the comparison between FEM and XFEM is fairly performed in terms of the approximation properties of the discrete functional spaces. The integration of the elementary matrices and vectors is therefore performed using the same quadrature. In that sense, the version of FEM used here is nonstandard because it also splits the multiphase elements for integration purposes in order to obtain the stiffness matrix and the force vectors with the same accuracy as the ones typically obtained with XFEM.

Example 1b (tilted interface):

We now consider the same problem as before, but with a slightly tilted interface (slope -15 • ). In that case, the solution is really 2D, and is plotted in Figure 3. As before, the FEM solution for a mesh that does not follow the interface is a rather poor approximation of the exact solution (computed using the FEM with a fitted mesh), both in terms of global energy norm and local evaluation of fluxes in the vicinity of the interface. Again, the XFEM behaves well in terms of global energy norm, practically coinciding with the exact energy, see [START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF] for analytical expressions of the same type of problems. However, this time, the XFEM seems to behave as bad as the FEM in terms of fluxes. This is due to the fact that the gradient of the solution may indeed be discontinuous, but that the chosen degrees of freedom, and the corresponding functional basis, do not allow to fix the gradients independently on each side of the interface. There is an additional freedom with respect to the FEM, but the exact solution is not necessarily contained in the approximating functional space. When the contrast in conductivity becomes large, large variations in the fluxes close to the interface can happen with small variations of the mesh.

Explicitly Enforcing Flux Continuity

We come back in this section to the way the flux continuity is enforced in the FEM and XFEM and propose to enforce it explicitly to improve the numerical results in terms of fluxes.

Weak approximation of Neumann Boundary Conditions and flux continuity on the interface

using (1c). In a similar way, integrating by parts in Ω 1 and Ω 2 and adding two expressions, a flux jump term appears in the integral form

Γ (q 2 -q 1 ) • n v dΓ (26)
that is readily taken as 0 attending to (1e). The numerical flux on Γ N converges asymptotically to g N as h tends to zero. The discrete approximation of the fluxes is such that equilibrium is reached in the discrete form, for the sum of the nodal versions of the body loads f and the surface loads g N . Thus, for a given mesh, the numerical flux on Γ N may be a poor approximation to g N , that is, it does not fit the data. The same occurs with the flux jump on Γ and, consequently, with the values of q 1 • n and q 2 • n which are, in practice, far from being accurate approximations.

Enforcing the flux continuity explicitly

In order to improve the quality of the fluxes on Γ, we propose to enforce explicitly the weak condition of continuity. That is, we propose to add the following condition to the XFEM solution, u X , b(µ, u X ) = 0, for all µ ∈ ṼH ,

where ṼH is a discrete space of weighting functions, that has to be chosen, and

b(µ, u) := Γ (q 2 -q 1 ) • n µ dΓ = Γ (ν 2 ∇u| Ω2 -ν 1 ∇u| Ω1 )•n µ dΓ. (28)
Note that the weighting functions µ should belong to L 2 (Γ). With this additional constraint, the flux continuity (1e) is expected to be fulfilled more accurately because it is not anymore competing with the residual coming from the interior of the domain.

The condition ( 27) is added to the original weak form (6) (or its discrete counterpart (10)) using the Lagrange multipliers approach. Thus, the resulting mixed problem reads: find u X ∈ V X and λ H ∈ ṼH such that

a(u X , w) + b(λ H , w) = ℓ(w) ∀w ∈ V X,0 (29a) b(µ, u X ) = 0 ∀µ ∈ ṼH (29b)
The matrix form of ( 29) results in an enlarged algebraic system of linear equations  

K uu K ua R u K T ua K aa R a R T u R T a 0     u a λ   =   F u F a 0   ( 30
)
where the matrices R u and R a are defined below. The Lagrange multipliers space ṼH must be chosen with care. If its dimension is too small, the flux continuity condition may be not be enforced strongly enough, so that the improvements might be limited with respect to the classical XFEM. If its dimension is too large, the resulting discrete problem may become unstable.

Indeed, we have observed through numerical tests that the obvious choice of selecting ṼH as the trace space (over the interface Γ) of the FE space V H does not introduce enough constraints to obtain a noticeable improving of the solution. Note that for this choice ṼH is generated by functions as the one depicted in Figure 4 and its dimension n I is equal to the number of element edges crossed by the interface.

An alternative choice, already considered in [START_REF] Soghrati | An interface-enriched generalized FEM for problems with discontinuous gradient fields[END_REF], consists in taken ṼH as the broken space generated by the functions Ñk depicted in Figure 4. The dimension of ṼH is 2n I because every edge crossed by the interface introduces two functions Ñk and Ñk+1 . For every k = 1, . . . , 2n I , the generic components of the matrices R u and R a in (30) read

[R u ] ik = b( Ñk , N i ) for i = 1, . . . , n H , (31) and
[R a ] jk = b( Ñk , RN j ) for j ∈ N a . (32)
We will show in Section 6 that the corresponding discrete problem is stable. Moreover, the flux jump on the interface is computed in a natural way in XFEM, leading to a simple computation of the Lagrange multipliers matrices R u and R a as

[R u ] ik = Γ (ν 2 ∇N i -ν 1 ∇N i ) •n Ñk dΓ and [R a ] jk = Γ (ν 2 ∇(RN j )| Ω2 -ν 1 ∇(RN j )| Ω1 ) •n Ñk dΓ
where the normal vector n is computed as defined in [START_REF] Chessa | An extended finite element method for two-phase fluids: flow simulation and modeling[END_REF]. As a final remark, and as is classically done in XFEM, it should be noted that the nodes are enriched only when the value of the level set in these nodes (φ i ) is, in absolute value, larger than a given threshold ǫ. Otherwise, the enrichment would only result in redundant parameterization. In all the computations performed in this paper, the threshold has been set to ǫ = 0.05, relative to a characteristic size of the smallest element of the mesh. On the other hand, it has been observed in the numerical experiments that normalizing functions Ñk in the broken space ṼH , for instance enforcing a unitary slope as in some implementations of the classical finite elements, does not bring a substantial improvement of the condition number of the global matrix.

Stability analysis

The method of Lagrange multipliers leads to a saddle-point problem. The infsup condition (or LBB) states [START_REF] Brezzi | On the existence, uniqueness and approximation of saddlepoint problems arising from lagrangian multipliers[END_REF] that the resulting formulation is stable if there exists a k > 0, independent of the mesh size, such that inf

µ∈ ṼH sup v∈V X b(µ, v) ||µ|| ṼH ||v|| V X ≥ k (33)
Note that the two norms in the denominator are different one from the other. One is defined over the domain Ω while the other is defined over the interface Γ.

It is often difficult to prove analytically that the LBB condition is fulfilled. In appendix A, we give one sketch of such proof of stability. However, as it is based on some strong hypotheses, that are difficult to prove in practice, we concentrate in this section on the Chapelle-Bathe test [START_REF] Chapelle | The inf-sup test[END_REF] that consists in solving the eigenvalue problem for a series of meshes and showing the boundedness of the smallest eigenvalue. Passing the numerical test does not guarantee stability but experience shows that it is very sensitive, and indeed manages to discriminate between stable and unstable problems.

The functions µ(x) ∈ ṼH and v(x) ∈ V X can be written as vector products µ(x) = Ñ(x) T µ and v = N(x) T v, where Ñ and N are vectors of the respective interpolation functions. The length (number of terms) of vector N is n x = n H + card N a = dim V X and the length of Ñ is 2n I . Their generic terms read:

[ Ñ(x)] i = Ñi (x), for 1 ≤ i ≤ 2n I , [N(x)] i = N i (x), for 1 ≤ i ≤ n H , and [N(x)] n H +i = R(x)N i (x), for i ∈ N a ,
and where µ and v are the vectors of the nodal values of the functions. The matrix form of the operator b(µ, v) in [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] can then be expressed as

b(µ, v) = Γint µ ν∇v • n dΓ = µ T Rv ( 34
)
where R = [R T u R T a ] is a 2n I ×n x rectangular matrix, whose blocks are defined in Equation [START_REF] Kim | A mortared finite element method for frictional contact on arbitrary interfaces[END_REF]. Assuming that, for every possible Lagrange multiplier µ ∈ ṼH , there exists a solution w ∈ V X , such that µ = ν∇w • n , we introduce the matrix J such that µ = Jw. This matrix J = [J T u J T a ] represents the flux jump operator in the bases Ñ and N. The general term [J] ij is such that where M µ is a mass matrix (corresponding to the L 2 product) for the elements of ṼH , with general term [M µ ] ij = Γint Ñi Ñj dΓ. We then obtain the discretized LBB condition:

2n I j=1 [J u ] ij Ñj = ν∇N i •n for 1 ≤ i ≤ n H , and 2n I j=1 [J a ] ij Ñj = ν∇(RN i )•n for i ∈ N a . We note that R = M µ J (35)
inf w∈R 2nx sup v∈R 2nx w T J T M µ Jv (w T J T M µ Jw) 1/2 (v T M v v) 1/2 ≥ k, ( 36
)
where M v is the mass matrix for V X . Note that all the norms are equivalent here, and we chose the L 2 norms for both || • || ṼH and || • || V X . The Chapelle-Bathe numerical test [START_REF] Chapelle | The inf-sup test[END_REF] consists in checking with numerical experiments the boundedness (from below) of the smallest non-zero eigenvalue of the generalized eigenvalue problem: find w ∈ R 2nx and λ ∈ R such that

J T M µ J w = λM v w. (37)
We therefore have to compute the smallest non-zero eigenvalue for a series of discretizations of the same domain with increasingly smaller elements. If that value does not appear to vanish for refined meshes, the test is considered passed. We consider for the test the Example 1b, described in Section 4 (see also Figure 6 for the solution obtained with our approach). The results are plotted in Figure 5 and indeed indicate that the test is passed.

Numerical examples

The examples of Section 4 show that, even in very simple scenarios, the FEM introduces a jump in the flux across the interface. Less intuitively, these examples show that the XFEM solution presents continuous fluxes over the interface only when the flux is orthogonal to the interface (see Figures 2 and3). In this section, we first show that our method works well for Example 1b. We then compare the behavior of our method to the FEM and XFEM solutions in the case of an unstructured grid. Finally, we describe the convergence of the solutions obtained with the three methods in terms of the norm of the flux along the interface. This convergence study is performed for a problem with analytical solution. For easier reference in the remainder of this section, we denote our method by XFEM+.

Example 1b.

We first apply the XFEM+ to Example 1b, which was already solved with the FEM and XFEM. The setting is described in Section 4 and maps of fluxes computed with the FEM and XFEM are presented in Figure 3. As illustrated in Figure 6(a), the fluxes computed with the XFEM+ are visually much more reasonable than those computed with the FEM or XFEM. Also, as expected, the flux jump cancels along the interface for the XFEM+. Note also that the flux jump is constant along the interface within one element with the FEM, while it is linear with the XFEM.

In Table 1, we present a comparison of the solution of Example 1b obtained with the three schemes. The comparison is performed both in terms of a global quantity, the energy norm, and local quantities, the maximum and mean flux jump along the interface. As explained earlier, the best solution in terms of the global energy norm is the XFEM, but the XFEM+ performs almost equally well. In terms of the local quantities of interest, however, the XFEM does not behave better than the FEM, and only the XFEM+ gives good results. In this example we analyze the value of the flux jump, rather than the actual value of the flux which could be analyzed using the reference solution for Example 1b provided in [START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF]. This is because we postpone this discussion to Section 7.3 in which the analytical solution is available and very simple to describe.

Unstructured meshes

In this example, we consider the same problem as before, with two small modifications: the slope is now -5 • and a non-structured mesh is used. The objective is to show that the results obtained earlier are not dependent on the specific alignment of the interface with the edges of the elements. The results obtained with the three schemes are presented in Figure 7. The fluxes in the elements crossed by the interface are plotted at the Gauss integration points, and the flux jump is plotted along the interface. The results confirm the conclusions of the previous example and imply that the method is valid whatever the relative alignment of the interface with the sides of the elements it crosses.

The proposed scheme adds restrictions to the XFEM functional space where the solution is sought. Consequently, the XFEM+ space is strictly contained inside the XFEM space. Thus, if measured in the global energy norm, the error for XFEM+ is necessarily larger than for XFEM. This is observed in the examples and shown in Table II. Nevertheless, the loss of accuracy associated with this phenomenon is very small (the energy norm of the XFEM+ error is only 0.4% larger than the XFEM one) if compared with the difference between FEM and XFEM (discrepancies of around 15%). In other words, a little price is paid in terms of energy to improve the accuracy of the fluxes.

Convergence study

The previous examples indicated that the flux jumps computed by the FEM and XFEM were not correct, and showed that the fluxes evaluated with the FEM and XFEM presented variabilities not compatible with the simplicity of the problems considered. However, we have not yet shown quantitative arguments to prove that the fluxes are indeed well estimated with the XFEM+. We therefore Let us consider the open disc Ω ∈ R 2 centered on (0, 0), and with radius 1. It is split into two concentric subdomains Ω 1 (disc) and Ω 2 (ring) at radius 1/2 (see Figure 8). We use cylindrical coordinates (r, θ), and define the parameter ν 1 = 1000 and ν 2 = 1. We finally consider the following problem: find u such that ν∆u + 4 = 0 , in Ω u = 0 , at r = 1 [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] with continuity of u and normal flux q = ν∇u • e r at the interface r = r 0 . The exact solution u ex of that problem is given by

u ex (r ≤ r 0 ) = r 2 /1000 -1 + 999/4000, u ex (r ≥ r 0 ) = r 2 -1. (39)
The exact value of the flux, equal in both Ω 1 and Ω 2 is

ν∇v = 2re r , (40)
so that the exact value of the normal flux along the interface is q ex = 2r 0 We compute the approximate solution to this problem using five different meshes with increasing number of elements and the three different schemes (FEM, XFEM, XFEM+). For each mesh and each scheme we compute the relative error in terms of the flux at the interface as

e = |q ex -q h | |q ex | , (41)
where q h is the approximate flux, computed on the interface Γ on the side of Ω 2 . We then plot in Figure 9 the evolution of that error with the size of the elements. This figure shows that the only method that converges with the size of the elements is the XFEM+. The convergence rate varies between about 0.5 and 1, although it is difficult to estimate the theoretical rate of convergence that the curve should follow. Interestingly, it can be noted that, in terms of the flux, the solution of the standard FEM deteriorates as the elements get smaller. Note, however, that we are measuring here only the accuracy of the flux close to the interface. Obviously the error in terms of the global energy norm gets smaller as the mesh is refined. The errors of the XFEM seem, at least for this example, to be independent of the mesh size and are maintained within the same order as h gets smaller.

CONCLUSIONS

In this paper, we have presented a modification of the XFEM for bi-phasic diffusion problems, that improves the evaluation of the fluxes in the vicinity of the interface. It requires the solution of a mixed problem and enforces weakly the continuity of the fluxes over the interface. The stability of the mixed problem has been discussed based on the Chapelle-Bathe numerical test, which was passed satisfactorily. For a very limited loss of accuracy in terms of global energy norm (in the worst case found in the examples, the undervaluation is of 0.5%, in the rest of cases, below the resolution of the four significant digits display), the examples presented in the paper showed that our method behaves extremely well in terms of the evaluation of the flux over the interface. The use of this improved method seems extremely appealing for the wide range of problems where the evolution of the interface is controlled by the fluxes along that interface. Note that the errors in the evaluation of the fluxes with the FEM and XFEM are particularly important when there is a large difference in properties between the two domains, so that our approach is particularly interesting in those cases (in the examples in this paper we have considered ratios of properties of 1000).

A Analytical study of stability

The numerical demonstration of stability presented in section 6, based on the ideas of [START_REF] Chapelle | The inf-sup test[END_REF] is here complemented by a theoretical result. This proposition states that for the proposed formulation, the LBB condition is proven to hold. This result requires some assumptions that are not obvious to guarantee but that, in practice, are fulfilled.

The proof is based on one idea already introduced in the numerical test. For the selected spaces, it is shown that any function µ in the Lagrange multipliers space, ṼH can be recovered as the jump of the normal flux of some function v in V X . Moreover, it is also shown that the norm of v is bounded by above by the norm of µ.

First, these properties are stated in an element-wise format.

Lemma A. = Γ e µ 2 dΓ. The standard FE shape function corresponding to the node P 1 is denoted N 1 , and the ridge function, defined at equation [START_REF] Dolbow | A numerical strategy for investigating the kinetic response of stimulus-responsive hydrogels[END_REF], is denoted R.

Then ∃α > 0 such that ∀µ ∈ Ṽ e H , ∃v ∈ span{N 1 , RN 1 } ⊂ V e X (i.e. describing v with the d.o.f. corresponding to P 1 only) verifying

1. ν∇v • n = µ; 2. v V e X ≤ α µ Ṽ e H .

Let us consider a function

v = u 1 N 1 + a 1 N 1 R ∈ span{N 1 , RN
:= ν∇N 1 R • n .
Note that γ 1 is a constant function on Γ e , which vanishes only if ∇N 1 is orthogonal to n, that is in the case of being Γ e perpendicular to the side P 2 P 3 of the element. In this case, the node P 1 has to be replaced by P 2 . Function γ 2 is linear (not constant) and therefore {γ 1 , γ 2 } is a basis of Ṽ e H . Thus, varying the coefficients u 1 and a 1 , the jump of v, may reproduce any function µ

∈ Ṽ e H , that is µ = u 1 γ 1 + a 1 γ 2 = ν∇v • n .
Moreover, being both v and µ represented by the coefficients [u 1 a 1 ] in the bases {N 1 , N 1 R} and {γ 1 , γ 2 } respectively, their norms are expressed as

v 2 V e X = [u 1 a 1 ]M v [u 1 a 1 ] T and µ 2 Ṽ e H = [u 1 a 1 ]M µ [u 1 a 1 ] T ,
being M v and M µ 2 × 2 mass matrices with generic terms

[M v ] ij = Ω e N 2 1 R i+j-2 dΩ and [M µ ] ij = Γ e γ i γ j dΓ for i, j = 1, 2
The second part of the statement is readily proved by taking α equal to αmax , the maximum eigenvalue of the following generalized eigenvalue problem.

M v q = αM µ q, (42)
The fact that this eigenvalue is bounded from below (it cannot be indefinitely small) is guaranteed by the fact that the elements are regular enough and the minimum length of the interface inside the element is also limited. This bounding property holds for a given value of H. It is necessary to check also the tendency of α in the limit case, that is when H tends to 0. Note that, with respect to the characteristic element size, H, M v scales with H 2 and M µ scales with H. Therefore, α is expected to scale with H. This is the right tendency for the final result shown in Theorem A. [START_REF] Vennat | Collagen fiber network infiltration: permeability and capillary infiltration[END_REF] because it has to be guaranteed that k is larger than a given value and k is taken as the inverse of α. Thus, taking k equal to the inverse of the α corresponding to H max (the largest possible element size) we are in the safe side.

The next proposition states that the decay of a function v defined in an element Ω e is fast enough to control the propagation of the norm into the neighbor elements.

Lemma A.2 (Bound of the norm propagated to the neighbor elements)

Let Ω e and Ω e+1 be two contiguous elements crossed by the interface and P 1 P 2 and P 3 the nodes of Ω e . Then, ∃β > 0 such that, for any v defined by the d.o.f. of Ω e , v ∈ span{N i , RN i } i = 1, 2, 3 [START_REF] Zlotnik | Numerical modelling of tectonic plates subduction using X-FEM[END_REF] it holds that

v V e+1 X ≤ β v V e X . (44)
Let us denote P 1 and P 3 the common nodes to Ω e and Ω e+1 , being P 2 the third node in Ω e , see figure 10. P 1 is selected such that it is on the same side of the interface as P 2 . The third node in Ω e+1 is denoted as P 4 as shown in Figure 10. In order to define a mapping transforming Ω e into Ω e+1 we introduce the following: let O be the intersection point of straight lines P 1 P 3 and P 2 P 4 , γ = |OP4| |OP2| and θ the angle between OP 1 and OP 2 . Then the affinity, A, of axis P 1 P 3 following the direction of P 2 P 4 and ratio -γ, maps Ω e into Ω e+1 . Selecting the proper coordinate axis, the Jacobian of the affinity is

J A = -γ -γ cot θ 0 1 ( 45
)
and the absolute value of its determinant yields

|J A | = γ. v 2 V e+1 X = Ω e+1 (u 1 N 1 + u 3 N 3 + a 1 RN 1 + a 3 RN 3 ) 2 dΩ = Ω e (u 1 N 1 + u 3 N 3 + a 1 RN 1 + a 3 RN 3 ) 2 Ω e+1 • A |J A |dΩ (46)
Note that due to the similarity of the shape functions N i Ω e+1 • A = N i Ω e for i = 1, 3.

Let us denote Γ as the inverse mapping of the interface Γ e+1 into Ω e , namely Γ = A -1 (Γ e+1 ). Note that the ridge function associated with Γ, R, is such that R = R Ω e+1 • A. The equation ( 46) can be rewritten as

v 2 V e+1 X = Ω e u 1 N 1 + u 3 N 3 + a 1 RN 1 + a 3 RN 3 2 γ dΩ. ( 47
)
Recall that The following lemma is a straightforward corollary of the two previous lemmas.

Lemma A.3 (Extension of lemma 1 to the entire space assuming the limited propagation of the Under the hypotheses of the previous lemmas, ∃α > 0 such that ∀µ ∈ ṼH , ∃v ∈ V X verifying 1. ν∇v • n = µ;

2. v V X ≤ α µ ṼH .
The following theorem is a direct consequence of the previous lemma and guarantees the fulfillment of the LBB condition.

||µ|| ṼH ||v|| V X = ||µ|| ṼH ||v|| V X
Since ||v|| V X ≤ α||µ||, taking k = 1/α the proposition follows. This is equivalent to satisfy the LBB condition and guarantees the stability of the scheme.

Support from Ministerio de Educación y Ciencia, Grant DPI2011-27778-C02-02 is deeply acknowledged.

Figure 1 :

 1 Figure 1: Illustration of (a) the complete domain Ω, the subdomains Ω 1 and Ω 2 and the interface Γ; (b) zoom on a generic point P of the interface Γ, with the definition of the normal and tangent vectors; and (c) description of the interface within an element.

Figure 2 :

 2 Figure 2: Example 1a (horizontal interface): the XFEM (right figure) performs better than the FEM (left figure) both in terms of global energy norm and local fluxes. The fluxes are evaluated at the Gauss points of each element and drawn in bold for the elements that are cut by the interface.

Figure 3 :

 3 Figure 3: Example 1b (tilted interface): the additional degrees of freedom in the XFEM (right figure) induce an improvement with respect to the FEM solution (left figure), in terms of the global energy norm, but not necessarily an improved evaluation of the fluxes close to the interface.

Figure 4 :

 4 Figure 4: Representation of typical elements of the trace space of V H (right figure) and of the broken space ṼH (left figure).

Figure 5 :

 5 Figure 5: Chapelle-Bathe test: evolution of the smallest eigenvalue as a function of the size of the mesh.

  (a) XFEM+ solution, u XFEM+ = 1Flux jump along the interface

Figure 6 :

 6 Figure 6: Example 1b: fluxes and energy norm of the solution computed with the XFEM+ (left figure, to be compared with Figure 3) and flux jump along the interface computed with the FEM (upper right figure), the XFEM (center right figure) and the XFEM+ (lower right figure).

  (a) FEM solution, u FEM = 1.924 (b) XFEM solution, u XFEM = 1.824 (c) XFEM+ solution, u XFEM+ = 1Flux jump along the interface

Figure 7 :

 7 Figure 7: Unstructured mesh example: fluxes and energy norm of the solution computed with the three schemes (figures (a), (b) and (c)) and flux jump along the interface computed with the FEM ((d), upper figure), the XFEM ((d) center figure) and the XFEM+ ((d) lower figure).

Ω 1 Ω 2 ΓFigure 8 :

 128 Figure 8: Discretized domain and interface for the convergence study

Figure 9 :

 9 Figure 9: Convergence of the relative error of the flux (not the flux jump) with the size of the elements for the FEM (dotted line), the XFEM (dashed line) and the XFEM+ (solid line).

(u 1 N 1 +Figure 10 :

 1110 Figure 10: Illustration of two contiguous elements Ω e and Ω e+1 for lemma A.2

Theorem A. 4

 4 Under the assumptions of regularity of the mesh and of local length of the interface, ∃k > 0 such thatinf µ∈ ṼH sup w∈V H b(µ, w) ||µ|| ṼH ||w|| V X ≥ kFor a given µ, take v as indicated in Lemma A.3. Then b(µ, w) = Γ µ 2 dΓ = ||µ|| 2 ṼH and b(µ, v)

Table 1 :

 1 Example 1b: comparison, in terms of global energy norm, maximum and minimum flux jumps (in absolute value) of the solutions obtained with the three schemes.

	scheme FEM XFEM XFEM+
	global energy norm 2.237	1.892	1.900
	max flux jump 66.3	37.06	1.84e-12
	mean flux jump 13.4	5.0	4.4e-13
	scheme FEM XFEM XFEM+
	global energy norm 1.924	1.824	1.824
	max flux jump	5.5	4.22	3.9e-12
	mean flux jump	2.4	0.5	1.4e-12

Table 2 :

 2 Unstructured mesh example: comparison, in terms of global energy norm, maximum and minimum flux jumps (in absolute value) of the solutions obtained with the three schemes.

  1 (Element-wise property)Let Ω e be one linear triangular element crossed by the interface Γ. The restriction of Γ to Ω e is denoted Γ e . The nodes of Ω e are denoted P 1 , P 2 and P 3 , choosing the order such that P 1 and P 2 are on the same side of the interface. As classically done in XFEM, we assume that ∃ǫ > 0 such that |Γ e | > ǫ. The restrictions of the functional spaces V X and ṼH to Ω e and Γ e are denoted V e X and Ṽ e H , with their respective L 2 norms v 2

	V e X	= Ω e v 2 dΩ and µ 2 Ṽ e H

PROBLEM STATEMENT Let us consider an open bounded domain Ω ⊂ R 2 , partitioned into two subdomains Ω 1 and Ω 2 with different physical characteristics (see figure 1 for notations). The boundary ∂Ω of the global domain Ω is divided into two parts on which two different types of boundary conditions (Dirichlet or Neumann) will be applied: ∂Ω = Γ N ∪ Γ D , with Γ N ∩ Γ D = ∅. The interface between the two subdomains Ω 1 and Ω 2 is denoted Γ = (∂Ω 1 ∪ ∂Ω 2 ) \ ∂Ω.We consider the following problem, which could for instance model the tem-

The interface condition (1e) is weakly enforced, in the same fashion as the Neumann B.C. (1c). Recall that, when integrating by parts the weighted residual of (1a), the term including the flux on Γ N ,Γ N ν∇u•nv ds is replaced by Γ N g N v ds