Kun Qiu

Siyuan Huang

Qiongwen Xu

Member, IEEE Jin Zhao

Member, IEEE Xin Wang
email: xinw@fudan.edu.cn

ParaCon: A Parallel Control Plane for Scaling Up Path Computation in SDN

Keywords: software-defined networks, distributed algorithms, performance management

The fundamental tasks of the control plane in Software Defined Networking (SDN) are to customize forwarding policies for the data plane and to provide global network view for applications. The logically centralized control plane design brings benefits in terms of network programmability and can largely ease network management. However, it also increases efficiency concerns. One practical control plane challenge is path computation, because it can require a significant amount of computation load if the network scale is large and the path requests from applications are frequent. In this paper, our goal is to build a high-performance control plane for path computation using multiple controllers. Previous works attempt to improve control plane efficiency by balancing only the load for data plane behavior between multiple controllers. Going beyond conventional wisdom, we designed ParaCon, a solution we propose to speed up the control plane by distributing the load of path computation. We also address the consistency and synchronization overhead challenges related to ParaCon design. To the best of our knowledge, ParaCon is the first attempt that utilizes node parallelism in SDN path computation. We evaluated ParaCon using both Mininet and real-world clusters. Our results show that the path computing time of ParaCon can achieve a speedup of 10x over Floyd (used in POX) and Dijkstra (used in ONOS) baseline implementations for networks with hundreds of nodes.

I. INTRODUCTION

R ECENTLY, Software-Defined Networking (SDN) be- came a hot topic in both academia and industry. SDN brings the benefit of network controllability by separating the control plane from the data plane. With OpenFlow [START_REF] Mckeown | Openflow: enabling innovation in campus networks[END_REF], a southbound interface specified by the ONF [2], a single centralized SDN controller can control the behavior of all the switches in a network. Despite its advantages, the single centralized structure may lead to scalability and performance challenges [START_REF] Hu | Scalability of control planes for software defined networks: Modeling and evaluation[END_REF]. The SDN controller can easily suffer from performance degradation resulting from a rapid increase of the network scale. Thus, such a design cannot scale up for realworld deployments, which is one of the reasons why current SDN deployments do not seem to take place in carrier-grade networks.

Therefore, improving the SDN controller performance and scalability is one of the major concerns in realizing SDN benefits. To address this challenge, academia and industry have proposed a diverse range of methods. One approach is to optimize the system architecture of existing single-server controllers for higher efficiency. To achieve better execution efficiency, some high-level language based controllers were redesigned as low-level language based controllers [START_REF] Gude | NOX: towards an operating system for networks[END_REF]. Another approach is offloading (partial or full) load from the controller to switches [START_REF] Yu | Scalable Flow-based Networking with DIFANE[END_REF], [START_REF] Curtis | DevoFlow: Scaling Flow Management for Highperformance Networks[END_REF]. Alternatively, utilizing multiple controllers in managing the SDN provides yet another choice [START_REF] Tootoonchian | Hyperflow: a distributed control plane for openflow[END_REF]- [START_REF] Panda | Scl: Simplifying distributed sdn control planes[END_REF] and the distributed control plane designs are becoming increasingly popular. For example, ONOS [START_REF] Berde | ONOS: Towards an Open, Distributed SDN OS[END_REF] claims to be the first carriergrade distributed SDN controller, and Onix [START_REF] Koponen | Onix: A Distributed Control Platform for Large-scale Production Networks[END_REF] is used in Google's B4 [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF] as their distributed control plane. Moreover, a new control plane design, which is based on a database approach was recently proposed [START_REF] Davie | A database approach to sdn control plane design[END_REF].

Existing multi-controller solutions aim to tackle the scalability issue in handling more switches, which may come at a prohibitive cost for maintaining and computing network information. More specifically, the SDN controllers periodically monitor the links, the forwarding tables and compute the network information such as the routing path and spanning tree of the global network. Among the variety of network information to be computed and processed in SDN controllers, the routing path is the most important one. Performance of mainstream SDN applications such as traffic engineering and content delivery networking [START_REF] Wichtlhuber | An sdn-based cdn/isp collaboration architecture for managing high-volume flows[END_REF] is tightly coupled with the performance of path computation. Path computation in a large-scale network indeed requires a significant amount of computation load. However, the path computation method in most of the existing SDN controllers is unable to meet the performance needs. As an example, ONOS [START_REF] Berde | ONOS: Towards an Open, Distributed SDN OS[END_REF], a widely used SDN controller, cannot achieve the intended throughput in a 205-switch topology due to its inefficiency in path computation. It is reported that the current ONOS design, though with multiple controller instances, only uses one single controller instance to compute paths, thus not fully utilizing the potential of parallelism in path computation [START_REF] Yamanaka | demo]scaling-up flow path computation for a large number of virtual sdn/nfv infrastructures[END_REF].

A straightforward approach to scaling path computation is to leverage the multiple controller instances in a distributed SDN control plane. For example, a possible solution is to implement parallel Dijkstra or Floyd algorithms [START_REF] Cormen | Introduction to Algorithms[END_REF] on multiple controllers. However, naïvely applying node parallelism to these algorithms within the context of the SDN control plane has suffered from two major challenges in practice.

First, the computation overhead is sensitive to topology changes. Even when there are only minor changes in the topology, offline algorithms need to start path computation on the updated topology from scratch. More specifically, allpairs path computation on the updated graph involves roughly O(n 3 ) time complexity [START_REF] Cormen | Introduction to Algorithms[END_REF]. Due to the potentially large size of SDN networks, frequent or even moderate network changes will pose a substantial amount of computation overhead. Second, parallel path computation needs to synchronize path information between involved controllers. Existing algorithms for parallel path computation are usually based on the Bulk Synchronous Parallel Computing (BSP) model [START_REF] Goudreau | Towards efficiency and portability: Programming with the bsp model[END_REF]. Fig. 1(a) shows that blocking synchronization in the path computation reduces the efficiency significantly. Moreover, synchronous path computation does not adapt 'online' property in the network either. With the increase in the number of involved SDN controllers, the synchronization overhead will result in significant performance degradation. Fig. 1. Synchronous and asynchronous path computation. In synchronous computation, the synchronization process separates path computation into several steps, and synchronization happened between the steps. In asynchronous computation, synchronization happened during the computation.

In this paper, we propose ParaCon, a distributed SDN control plane that addresses the performance challenges of distributed path computation of SDN. Under this design, SDN programmers can settle path computation without worrying about the details of synchronization while achieving the benefits offered by parallel computing. To the best of our knowledge, ParaCon is the first control plane architecture to utilize parallelism in path computation for SDN. More specifically, ParaCon introduces several new features to make path computation in multiple controllers efficient. It introduces a new distributed model in SDN that is capable of asynchronous path computation like Fig. 1(b) between multiple controllers.

ParaCon includes an online algorithm that enables all SDN controllers to make incremental changes to the computed path results upon receiving a particular topology change while existing approaches need to restart the computation from scratch. It makes ParaCon more efficient in scenarios with frequent topology updates. Moreover, it seeks to minimize the overhead between multiple controllers by designing an asynchronous algorithm. In this algorithm, different controllers can synchronize information without blocking. Finally, it introduces a hybrid consistency model to maintain the topology-related information and path information between controllers during the asynchronous path computation. It utilizes both strong consistency and eventual consistency in the maintenance of this information for further improving the performance. Our theoretical analysis shows that the online algorithm has a bounded convergence time. Moreover, experimental results show that ParaCon can significantly reduce the path computation time when the scale of the network is increasing.

The rest of the paper is organized as follows. Section II gives an abstraction of the problem space and challenges. Section III gives an overview of ParaCon architecture. Section IV gives the details of the algorithm. We present our evaluation in Section V and related work in Section VI. We conclude in Section VII.

II. PROBLEM STATEMENT

As we have mentioned above, we seek to minimize the overhead between multiple controllers by designing an asynchronous path computation algorithm. To motivate the asynchronous computing approach, we need to consider several aspects, among which the consistency problem is the most important one. Namely, we need to consider what information is to be synchronized, when to synchronize it, and how to minimize the synchronization overhead during the path computation.

Generally speaking, to compute path information in a distributed fashion, we need to maintain at least the following information:

1) Link view: The topology-related information, i.e., the link/switch status. 2) Association information: The association between switches and controllers. 3) Path view: The path information in each controller during the path computation.

From the example described in Fig. 2, we can see the link view (topology-related information) needs to be synchronized between all controllers to make sure the applications in different controllers get the same global topology. Otherwise, an inconsistent link/switch status will make an inconsistent result of path computation. We refer to this synchronization policy as strong consistency policy.

Before we define the consistency policy formally, we introduce some notations. Let H be some new information updated to a controller such as a link state or path change. For a given controller c ∈ C, where C is controller set, p c indicates its status such as p c ∈ {N O, Y ES}: when in status N O, a controller has not received or accepted H, when in Y ES it has received or accepted H. The control plane is defined as in a consistent state after H only if ∀c ∈ C, p c = Y ES. The multicontroller SDN control plane is said to be available when it can reply to control plane request messages. Otherwise, it is said to be not available. Different consistency models impose diverse levels of strictness on the availability of the control plane during in the inconsistent state. A strong consistency policy only allows new changes to be accessed after they have been applied to all controllers. Fig. 2. Suppose there are four switches S1, S2, S3 and S4 in a network. Controllers C1, C2 are added to the network with the association. The topology-related information needs to be synchronized to controller C1 and controller C2 to make sure they have the same path view on switches S1, S2, S3 and S4. Otherwise, an inconsistent link/switch status will make an inconsistent result of path computation. For example, if we do not synchronize the link status when the link between S1 and S2 fails (we use a red line between S1 and S2), we may get an inconsistent path view (the path from

S1 to S2 is S1 → S2 → S3 in C1, but S1 → S3 in C2).
Meanwhile, association information also needs to be synchronized to each controller with a strong consistency policy, or the assignment between controllers and switches may be incoherent if migrations of switches happened in the meantime.

However, synchronizing path view (path information) during the path computation, even at a minimum required frequency, may eventually lead to network congestion due to the excessive overhead in control message exchange under a strong consistency policy [START_REF] Decandia | Dynamo: amazon's highly available key-value store[END_REF]. For example, Two-Phase Commit [START_REF] Mohan | Transaction management in the r* distributed database management system[END_REF] or Paxos [START_REF] Lamport | Paxos made simple[END_REF] are revealed to be inefficient by signaling overhead.

Therefore, it is necessary to reduce the synchronization overhead of path information. Considering that the three types of information we mentioned above involve different transmission overhead, we need to study the amount of exchanged messages quantitatively. We assume the exchange process on the control plane adopts a message-passing model [START_REF] Canini | A distributed and robust sdn control plane for transactional network updates[END_REF] [23], and we consider the minimum number of messages transmitted for reflecting the change in a link/switch status and a switchcontroller assignment is 1 unit. We have two types of operations: synchronizing topology-related information and synchronizing path information. Synchronizing topology-related information involves one simple message that is associated with one switch. In contrast, synchronizing path information requires more transmission overhead. If there are n switches, the controller needs to transmit up to n 2 messages because the change is associated with the path information that has n 2 paths at most.

Accordingly, synchronizing the path information can be much more resource-intensive than other operations. This motivates us to minimize the overall overhead in synchronization by reducing message transmission for path information. Different from strong consistency policy, an eventual consistency policy is a specific form of weak consistency. The SDN control plane guarantees that if there are no new updates such as link status change, eventually all controllers will return the same path information, which is computed based on the topologyrelated status. Consequently, an eventual consistency policy can be used for synchronizing path information to reduce the transmission overhead. Instead of synchronizing path information to all controllers immediately, an eventual consistency policy allows a delay in synchronization. The path information can be accessed during the time with inconsistent results. Naturally, the eventual consistency results in inconsistent path information during the inconsistency window. To avoid inconsistency effects as much as possible, each controller should independently compute the path information based on topology-related information to ensure consistency. However, with the increasing scale of the network, the computation overhead will inflate rapidly since the time complexity of allpairs path computation algorithm is O(n 3 ) in general [START_REF] Cormen | Introduction to Algorithms[END_REF].

In summary, finding an equilibrium between synchronization overhead and computation time is difficult. It motivates us to propose a distributed architecture, which neither significantly increases synchronization overhead nor obviously increases computation time for parallel path computation.

III. THE PARACON ARCHITECTURE

The core objective of ParaCon is to offer a highperformance SDN control plane for path computation with minimal transmission overhead between multiple controllers. Fig. 3 shows the architecture of ParaCon. Our architecture is based on a distributed structure with an asynchronous parallel algorithm. We will give a brief introduction to its two components.

A. The Distributed Structure 1) Switch: The data plane consists of switches and links. If a switch needs to forward a packet to other switches, it will request a path between two switches from the controller.

2) Topology Abstraction: We abstract the topology using a graph [START_REF] Raghavendra | Dynamic graph query primitives for sdn-based cloud network management[END_REF]. The global topology view includes both a link view and a path view. Link view is composed of nodes and edges, where nodes are switches and edges are links. Link view indeed reflects the topology of the underlying network. Path view provides information on the best route between any two switches. Path view indeed reflects a snapshot of the weights for the all-pairs shortest paths at a given time. As the algorithm progresses, the path view will be updated till the results at all the controllers converge.

3) Controller and Allocator: The control plane is distributed. Each controller monitors the status of their switches/links, which are allocated by the Allocator. If the status of switches/links has changed, the controller will synchronize the topology-related information to other controllers by utilizing the strong consistency policy. Thus, all controllers are sharing the same global topology view based on switches/links status. Meanwhile, each controller also computes the path information of their switches and receives the topologyrelated/path information from other controllers.

B. The asynchronous parallel algorithm

An asynchronous parallel algorithm is introduced, which uses the eventual consistency policy to improve the efficiency of the distributed controller design and eliminates the effects of the inconsistency due to the eventual consistency policy. More specifically, the design choice we choose is to synchronize the topology-related information (link view) and the affiliation with a strong consistency policy and to synchronize path information (path view) with an eventual consistency policy, with the goal of making a correct balance between transmitting overhead and computing overhead.

Furthermore, different from existing algorithms such as Floyd-Warshall [START_REF] Cormen | Introduction to Algorithms[END_REF] or Dijkstra [START_REF] Cormen | Introduction to Algorithms[END_REF], which are used in existing controllers, the proposed algorithm utilizes the 'online' property to reduce synchronization overhead. Upon a new topology change event, it does not need to wait until the end of updating the previous event. Instead, our algorithm can handle the new event immediately.

IV. PARALLEL PATH COMPUTATION

As mentioned above, we synchronize the link/switch status and the affiliation of switches with a strong consistency policy. To this end, we use Two-Phase Commit policy to submit the change so that all controllers can receive it. This section details the asynchronous path computation algorithm. The algorithm is more complicated than the synchronous one because of its interweaving with the eventual consistency policy. We use several steps to describe our algorithm.

A. Key notions

Switch and controller set: There are switches and controllers in a Software-Defined Network. The controllers are spatially distributed, and each switch is only controlled by one controller. Set S stands for the switches and C stands for the controllers. There are n switches s 1 , s 2 , s 3 , ..., s n in set S, and m controllers c 1 , c 2 , c 3 , ..., c m in set C. The element in vector Con(s) means the controller which controls switch s.

Link view: Link view is a graph L(S, E) where S stands for the switches set and E stands for the links set. All elements in set E are two-tuples (x, y), x, y ∈ S. The link status between two switches is expressed in an integer as the weight w. The link view is thus represented as an adjacency matrix. L a,b stands for the weight of the link between s a and s b . If there is no link between s a and s b , the L a,b = -1.

Path view: The path view provides information on the cost of the best (usually the shortest) route between any two switches. A path view can also be regarded as a weighted full graph T (S, W ), where S stands for switches and W stands for the weight of the shortest path between two switches. The weight of the shortest path is expressed as an integer, which is calculated based on link view L(S, E). The path view can be implemented as an adjacency matrix as well. T a,b stands for the weight of the shortest route between s a and s b . T a, * stands for the weights for all the shortest routes beginning with s a .

Forwarding table: Each switch has a list of forwarding rules to indicate the next hop of the best path to other switches. When we compute path view, we can also get the corresponding forwarding table for the shortest paths. The forwarding table is described as a matrix. F a,b represents the next hop from s a to s b , which is the first constituent edge along the shortest path from s a to s b .

In summary, we use 1) link view L as the topologyrelated information; 2) vector Con(S) as switch-controller association; and 3) path view T and forwarding table F as path information.

B. Centralized approach

We discuss a centralized algorithm first. There are several centralized algorithms for computing the paths in a network. The well-known Floyd-Warshall [START_REF] Cormen | Introduction to Algorithms[END_REF] and Dijkstra [START_REF] Cormen | Introduction to Algorithms[END_REF] algorithms have an O(n 3 ) time complexity for the all-pairs shortest path (APSP) problem in a network with n nodes. These algorithms have been used in various SDN controllers. For example, POX uses the Floyd-Warshall algorithm. However, both algorithms are offline algorithms, i.e. they can not deal with any link change while the algorithm is in progress. In contrast with existing approaches, we first design an online APSP algorithm based on Moore [START_REF] Cormen | Introduction to Algorithms[END_REF] algorithm and then tailor it for parallel scenarios. Note that Moore algorithm is an optimized version of the Bellman-Ford algorithm used for solving the single-pair shortest path (SPSP) problem.

The centralized algorithm uses a queue at each controller to store the nodes that will be checked for path cost update. We briefly explain how it works. If the queue is not empty, the algorithm will fetch the head of the queue and relax it. Relaxing a node involves updating the cost of existing paths if the total cost can be reduced by passing through the node.

If the node can be relaxed, which means that there is at least one path which becomes shorter by passing through this node, its neighbors will be inserted into the queue.

The correctness of our algorithm can be proved by contradiction. We can prove all-pairs shortest path are computed when the queue is empty. We use a contradiction to give simple descriptive proof: First, we assert the following Lemma: shortest path is composed of shortest path. In other words, for a shortest path x → . . . → a → . . . → b → . . . → y, any sub-path a → . . . → b must be a shortest path. The proof is based on the notion if there was a path shorter than any sub-path a → . . . → b, we can replace the subpath with the shorter path to make the whole path shorter. We want to prove the Proposition: all-pairs shortest paths are computed when the queue is empty. Suppose for the sake of contradiction that is not true: there exists one path that is not the shortest one when the queue is empty. We assume the path is x → . . . → u → . . . → y. With the Lemma we proposed above, we can find u that satisfy the part of the path x → . . . → u is the shortest path, and the other part of the path u → . . . → y is not the shortest path. Thus, we can always find another path u → v → . . . → y passing through v to make it shorter. To relax v, the algorithm will put u and its neighbor v into the queue, which is a contradiction to the Proposition.

Example: Suppose we have 4 switches and 1 controller. The link view and path view are shown in Fig. 4. When we change the link status between s 0 and s 1 to 30. The changes of topology view and the queue P and Q are shown in Fig. 4.

The approach can deal with link status change by just putting switches into the queue and waiting for the result. Because of its 'online' property, the time complexity is O(p * e * n), where e is the number of links and p is related to the diameter of the graph. In general, if the network diameter is small, e.g., the diameter of some topology in data-center is less than 3, the iterations will be completed very fast. In contrast, some extreme topologies with a large diameter made the p become large due to an iteration oscillation. For example, a topology with a large linear subgraph (in which most nodes have a degree ≤2) will lead to more iterations. Fortunately, the graph of real network topology rarely has a linear subgraph, and the oscillation can be with further optimizations.

The algorithm runs in a centralized way in each controller, which has no additional node-parallelism advantages compared to the other existing algorithms.

C. Asynchronous parallel path computation

Adopting a parallel path computation approach makes it all possible for all controllers to compute the path collaboratively. The right matrix shows the weight of shortest routes between switches. We use the blue square with the blue background to indicate that the path view is updated by the controller. After link 0, 2 has changed to 30, we put 0 and 2 to the queue Q. The algorithm gets the head of Q and computes the shortest route of it. After processing 0, we add the adjacent: 1 and 3 to Q because there are some changes when processing 0. If the queue is empty, the algorithm stops. In this example, the algorithm uses 6 steps: 0, 2, 1, 3, 0, 2 to finish the algorithm.

We assign the switches and computation load to different controllers. However, due to the heterogeneity in computing power and transmission overhead, the inconsistency of path view becomes severe under this design. A straightforward approach to consistency is to add a barrier to the queue after an iteration has been done. The barrier is used for synchronizing the path computation process for all controllers. If the head of the queue is a barrier, the path computation thread will be blocked until other controllers have reached the barrier. However, adopting a barrier will increase the synchronization delay substantially.

To turn our algorithm into an asynchronous one, we redesign the barrier to trade consistency for synchronization overhead. Thus, we employ the eventual consistency. A weak synchronization signal, or SYNC, is introduced which allows a controller to send its path view to other controllers without blocking any computation thread.

Another concern is that we need to deal with the rising conflicts when synchronizing the path view in different controllers with the eventual consistency. It is clear that some conflicts must be overwritten, and some conflicts cannot be overwritten. Otherwise, the algorithm may take extra time to compute the path (e.g., if a new path is overwritten by an old one), and it can significantly increase the overhead in our control plane. We introduce the following policy for path overwriting.

1) If the received path information is newer than the existing one, the existing one needs to be overwritten 2) If the received path information is older than the existing one, the existing one need not be overwritten The chronological order, either newer or older, is based on the logical clock. This motivates us to use Vector Clock to establish the global logical clock for each element in path view T .

Besides newer conflicts and older conflicts, there still exist other conflicts, which are neither newer nor older in vector clock that should be solved manually. Fortunately, under ParaCon design, these conflicts can be overwritten immediately, because directly overwriting path view is equivalent to changing the sequence of the queue Q in the centralized approach. This operation does not affect the results of the path computation which is an iteratively based approach.

Algorithm 1 is the ParaCon parallel path computation algorithm we proposed; it keeps the path view eventually consistent. The algorithm is deployed into all controllers and makes controllers compute the paths concurrently. The algorithm has a loop and is separated into three parts. First, it fetches the head of the queue and checks if the head is a synchronization signal. If yes, it sends view other controllers and continues the loop. Otherwise, it relaxes its path view of the head and puts the adjacency nodes which are relaxed by the head to the queue. At last, it pushes a synchronization signal into the queue. The algorithm can do path computation in multiple controllers efficiently. During the inconsistency window, the controller can fall back to ondemand mode (utilizing Dijkstra/Floyd algorithm) to obtain correct results.

Example: Suppose we have a topology with 4 switches and 2 controllers. Switch s 0 and s 1 are controlled by controller c 0 , while s 2 and s 3 are controlled by c 1 . When we change the link status (weight) between s 0 and s 2 to 30, the changes of path view and the queue Q 0 , Q 1 are shown in Fig. 5.

D. Convergence analysis

We use an epidemiology model [START_REF] Jelasity | Epidemic-style proactive aggregation in large overlay networks[END_REF]- [START_REF] Holliday | Epidemic algorithms for replicated databases[END_REF] to analyze our algorithm. Usually, the epidemiology model is used as an analyzing tool for network protocols such as Gossip or RIP (Routing Information Protocol) [START_REF] Boyd | Randomized gossip algorithms[END_REF], [START_REF] Haridasan | Gossip-based distribution estimation in peer-to-peer networks[END_REF].

Firstly, we define the symbols. We suppose there are n nodes (switches) in the graph (link view), and we assume k as the average number of nodes relaxed in the graph when SYNC is received (from line 5 to line 9 in Algorithm 1 ). To put this another way, when SYNC is received, the cost of existing paths can be updated as they can be shorter through these k nodes. Thus, we define β as update rate with β = k n , (0 < β < 1). As we have mentioned in the centralized approach, each node has two states: it can be relaxed or cannot be relaxed. Nodes will stay in the second state (cannot be relaxed) eventually and will not convert into the first state (can be relaxed). Fig. 6. An example for explaining the convergence analysis. There are n = 9 nodes in the graph. Each node has two states: it can be relaxed (red) or cannot be relaxed (green). There are 4 nodes relaxed at time t 1 , but only 3 nodes stay in the second state (the left one may have chances to be relaxed by other nodes). At time t 2 , 2 nodes relaxed in the graph, and both of them stay in the second state. Finally, all nodes will be in the second state eventually at the last time te. We assume kn as the number of nodes relaxed at time tn. In this example, at least we have k 1 = 3, k 2 = 2. k is the average number of nodes relaxed in the graph, which is the average number of k 1 , k 2 , . . . , ke.

Obviously, k is closely related to the topology, such as the average degree of nodes or the diameter of the graph.

At any time t, we assume the number of nodes in the first state is y t and the number of nodes in the second state is x t . Thus, x 0 = n, y 0 = 1 when a link status change happened, and at any time, x t + y t = n + 1. Also, we have dx dt = -βxy and we get

y = n + 1 1 + ne -β(n+1)t
Obviously, lim

t→∞ y = n + 1
Thus, all nodes will be in the second state finally. We assume t = clog(n), we can get

y ≈ (n + 1) - 1 n ck-2
We give an example to describe these symbols in Fig. 6. As the update rate β is related to k, which is influenced by the property of the graph, i.e., the diameter and the average degree of nodes, the number of nodes that are not in the second state at a given time t, is widely ranged. For example, k in a full mesh topology is larger than k in a linear topology, and obviously, the convergence time of a graph with a small diameter is less than the time of a graph with a large diameter. We have confirmed this conclusion by evaluation results in section V.

V. THE PARACON PROTOTYPE

Our work is motivated by the intuition that the performance of multi-controller solution can be improved using parallel computing in a cluster. We believe the most exciting results from the efficiency of the structure can be achieved in a realworld implementation.

A. Implementation

We implemented ParaCon structure using commodity hardware. The distributed control plane is based on 10 Dell Pow-erEdge 2950 (Xeon E5405 with 2G RAM) servers. We have four virtual machines in one server, and each virtual machine handles one controller. Therefore, we have 40 controllers at most. The controller is based on a modified POX, in which we replace the module of path computation by our parallel version. The switches are provided by Mininet [START_REF] Lantz | A network in a laptop: rapid prototyping for software-defined networks[END_REF]. We convert topologies into a Python program, which is used as topology input by Mininet. Due to its programmability, we can use any topology to test our distributed control plane. We will discuss the details as follows.

1) Path View Module: The path view module is separated into three threads: computing the path, handling event of topology-related information change, transmitting and receiving messages. The module uses an adjacency matrix to store the path view. We also implement a Two-Phase Commit to submit the change of link/switch status and assignment of switches. The transmitting and receiving thread uses TCP to communicate with its counterpart in other controllers. To improve efficiency, we employ eventlet, a concurrent networking Python library that provides highly scalable non-blocking I/O [START_REF] Eventlet | [END_REF].

2) Modified POX: We implement ParaCon based on POX, a widely used controller in SDN research community. Two modules in POX are modified, namely, Openflow.Discover, which controls the topology discovery, and Forwarding.l2_multi which controls path computing. Briefly, first, POX catches link change event by LLDP protocol using Openflow.Discover, and send to path view module. Next, our module sends an event of topologyrelated information change to all controllers and computes the path using parallel computing. After the computing is finished, we get the route and modify the flow table in OpenFlow switches using module Forwarding.l2_multi.

B. Processing path computation after link change

To evaluate ParaCon's performance in processing link change, we generate several full mesh topologies ranging the number of switches from 50 to 300. With the dense topologies, we can balance the heavy path computation load to controllers using even switch assignment. Meanwhile, we also change the number of controllers from 1 to 40 and assign switches to controllers. Even if the common number of controllers is generally less than 10 [START_REF] Berde | ONOS: Towards an Open, Distributed SDN OS[END_REF], we still set up to 40 controllers in order to test its scalability. We assign each controller 1 CPU core. We randomly shut down links and measure the time from when the link fails to when the new all-pairs path computation is over. In Fig. 7 we can see that by increasing the number of controllers, the computing time is reduced significantly. However, with 40 controllers, each (PowerEdge 2950) server has to run 4 controllers, each occupying 1 CPU cores; since each server only has 4 CPU cores, and the hypervisor uses at least 1 core for scheduling, the performance of each controller gets decreased. This explains the slight improvement in computing time between 20 controllers and 40 controllers. We also compare our module to other existing modules: a standard POX module (implemented by Floyd-Warshall algorithm), and an OSPF module (implemented by a distributed Dijkstra algorithm). From Fig. 8, we see that with an increasing number of switches, the computing time used by ParaCon is far less than other modules.

We count the average frequency of communication between controllers in one controller. From TABLE I, we can see that with an increasing number of controllers, the communication overhead is well balanced into all controllers.

C. Comparing with ONOS

1) Setup: To study the performance gain of ParaCon over the state-of-the-art solution, we compare it with ONOS v1.3, a popular and stable distributed OpenFlow controller [START_REF][END_REF]. Note that the path computation in ParaCon is in a different way than is done in ONOS. ONOS uses an ondemand mechanism of path computation, but ParaCon uses a pre-computed one. In ONOS, the forwarding module (fwd.ReactiveForwarding) receives the PacketIn and fetches the source and destination address. Then it calls the path computation modules to get a path from the source to the destination. ONOS provides a variety of algorithms for path computation. The default algorithm getPaths() uses a DijkstraGraphSearch() method. Finally, the forwarding module converts the path to flows and inserts them into switches. Therefore, we compare the path computation module of ParaCon with the getPaths() (using Dijkstra by default) method in ONOS.

In this evaluation, we do not use full mesh topology because a Dijkstra-based module (e.g., OSPF, ONOS) is obviously slower than ParaCon with an increasing number of nodes in a full mesh topology. Therefore, we compare path computation module between ParaCon and ONOS using real topologies. We set up 4 ParaCon controllers, and we use three typical topologies as shown in TABLE II. 'CERNET' and 'USCarrier' are provided by Topology Zoo [START_REF]Topology zoo[END_REF]. In order to use these topologies, we first get topology files from the Topology Zoo. These files are on a 'GML'/'GraphML' format, and we use Python and its library iGraph to convert these files into an adjacency matrix and edge list. For the topologies we made by ourselves ('Leaf-Spine'), we directly generate the adjacency matrix and edge list. 'CERNET' is a small topology with a small diameter; while 'USCarrier' is a large topology with a larger diameter. The third one 'Leaf-Spine' is a generated topology, which is usually used in a data center. Although it has a large number of nodes and links, the diameter is the smallest (only two hops from any hosts to any hosts in this topology). 2) Path Computation Latency: It is not straightforward to evaluate the performance between ParaCon and ONOS due to the different mechanisms for path computation (pre-compute v.s. on-demand). If there's no topology change, ParaCon can reply to the path request immediately with no extra latency because the results have been pre-computed. Consequently, ParaCon will outperform ONOS in scenarios with infrequent topology change.

For a fair comparison, we consider the case with link failures. Suppose the network topology has changed as a result of a link failure, the controller must re-compute the path. We compare the average latency for answering a batch of path requests. For ParaCon, the latency mainly relies on APSP computation despite the batch size, because ParaCon only computes APSP once, and it does not consume any time in querying a path. Therefore, to evaluate the computing time in ParaCon, we set up the topology and make a link fail, then we measure the running time from when the link fails to when new APSP are computed. Finally, we divide the running time with the number of requests to get the average computing time of one request. For ONOS, we set up the same topology and make a link fail; then we also measure the running time of getPaths() method for getting a path in the new topology. The running time for one request is stable, no matter how many requests it received. It is because the running time of getPaths() method only depends on the number of nodes (switches). Note that Dijkstra algorithm can have a time complexity of O(e + n log n) if implemented with a Fibonacci heap. However, since the original ONOS only utilizes the Dijkstra algorithm without Fibonacci heap, and the Fibonacci heap can take a very long time to complete some operations in the worst case [START_REF] Fredman | The pairing heap: A new form of self-adjusting heap[END_REF], we consider the general case O(n 2 ) [START_REF] Cormen | Introduction to Algorithms[END_REF].

From Fig. 10, we can see that ParaCon performs better than ONOS in 'CERNET' and 'Leaf-Spine' because these topologies have a small diameter. Although the iterations of our algorithm are closely related to the diameter of the graph, the increasing number of nodes and links does not reduce the performance of the algorithm significantly. The 'USCarrier' is the worst case for ParaCon: if the number of queries is smaller than 60, the efficiency of ParaCon is less than ONOS. However, with the increase in path requests, the computing time decreases substantially. In 'CERNET', the efficiency of ParaCon exceeds that of ONOS in the early stage. In 'Leaf-Spine', the ParaCon uses less time to finish the computation of APSP than ONOS, while ONOS uses more time to compute only a single-pair path. We believe ParaCon is more suitable for a data center environment than ONOS, given that data center topologies usually have a small diameter.

3) Recovery Time Upon Failure: We define recovery time as the time between when a link fails and when a new path is obtained. From TABLE IV, we can see that in 'CERNET' and 'USCarrier', ParaCon only increases the overhead slightly, but in 'Leaf-Spine', ParaCon is much faster than ONOS. This is because a larger number of iterations will make ParaCon perform worse. A topology with a large diameter will make the number of iterations larger. In the linear part of the topology, where nodes have fewer degrees than in the non-linear part of the topology, path change information propagates slowly. It will increase a larger number of iterations in ParaCon. The observation provides evidence for the result that there are many linear topology subgraphs in 'USCarrier'.

Although the new path upon link failure is computed, queues in controllers may not be empty. The controllers are still computing until queues are empty. In other words, at the time "Recovery Time Upon Failure", the new path upon link failure is computed, but not all the all-pairs shortest paths are fully computed. Thus, the recovery time is smaller than the computing time for the all-pairs shortest paths, as mentioned in the previous section.

4) Convergence Upon Incoming Requests from Switches:

As the eventual consistency application is one of our contributions, whose impact we want to assess precisely, we export the intermediate results during path computation to measure the 'dynamical correctness ratio'. The dynamical correctness ratio is defined as a percentage of consistent results, to show the consistency influence on the results during the computation. A lower correctness ratio may increase the possibility of getting a stale path. We used the same topologies ('USCarrier', 'CERNET' and 'Leaf-Spine') to perform this evaluation. For each topology, we randomly change the topology 30 times by selecting two switches, and randomly adding/removing the link between them with a 1s time interval. We measure the correctness ratio during the process. From Fig. 10, we can see that the results are confirming our intuition that a graph with small diameter ('Leaf-Spine') gets a smaller convergence time, and a graph with a large diameter (e.g., 'USCarrier' has a linear part) makes the convergence time larger. Comparing Fig. 9 (a) with Fig. 9 (b), we can see that most topology changes in 'USCarrier' make the correctness ratio go down to 80% ∼ 90%, and it uses more time to recover. But in 'CERNET', topology changes only make the correctness ratio going down to 95%, with a faster recovery than in 'USCarrier'. In 'Leaf-Spine', the correctness seldom changes, which is pretty stable around 100%, and the recovering time is the fastest. In summary, a graph with a large diameter (especially with some linear part) will lead to a large number of iterations, which increases the amount of data that needs to be synchronized, besides the convergence time.

5) Synchronization Data Overhead: One way to evaluate the impact of the communication overhead is to quantify the amount of data that needs to be synchronized. As defined in Section II, we consider the minimum number of messages transmitted (such as a link/switch status message) as 1 unit. As for the previous evaluation, we measure the amount of data synchronized in controllers with 30 random topology changes VI. RELATED WORK 1) Large-scale graph computation: As the scale of many graphs, e.g. web graphs or social network graphs, can reach billions of nodes and trillions of edges, data processing over such graphs becomes a concern. Some research such as GraphChi [START_REF] Kyrola | Graphchi: Large-scale graph computation on just a PC[END_REF] and BBQ [START_REF] Xu | Fast shortest-path queries on large-scale graphs[END_REF] enables large-scale graph computation using just one PC. Obviously one single PC has only limited scalability in large-scale graph computation. Thus, Pregel [START_REF] Malewicz | Pregel: a system for large-scale graph processing[END_REF] presents a parallel computation model for solving this task. Other solutions such as GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF] and PowerGraph [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF] also focus on increasing the performance of large-scale graph computation. However, most of them are BSP-based algorithms, which are inflexible in networking path computation.

2) Distributed shortest path algorithm: The shortest path block is the key one of many network applications, such as related to routing or failure recovery. Most of distributed shortest path algorithms are variations of either the Dijkstra or Bellman-Ford algorithm. J. K. Antonio [START_REF] Antonio | A fast distributed shortest path algorithm for a class of hierarchically clustered data networks[END_REF] proposed an O((n/p) log n) (n is the number of nodes, p is the number of processors for each node) algorithm for solving all-pairs shortest path problem on balanced hierarchically clustered (BHC) topology. Hence, S. Zhu and G. M. Huang [START_REF] Zhu | A new parallel and distributed shortest path algorithm for hierarchically clustered data networks[END_REF] proposed to adopt a newer algorithm for solving the multiple origins shortest path problems (MOSP) on hierarchically clustered topology, where m is an integer related to the topology. It is worth mentioning that there is also a work [START_REF] Chakaravarthy | Scalable single source shortest path algorithms for massively parallel systems[END_REF] designed for computing Single Pair Shortest Path (SPSP) in a High Performance Computer (HPC) with 500, 000+ CPU cores system. Moreover, some approximation algorithms can speed up the shortest path computation. Suppose D is the diameter of the graph and n is the number of nodes in the graph, B. Patt-Shamir and C. Lenzen [START_REF] Nanongkai | Distributed approximation algorithms for weighted shortest paths[END_REF] proposed a distributed approximation algorithm with Õ(n (1/2)+(1/2k) + D)-time, (8k log(k + 1) -1)-approximation (k ≥ 1) for single-pair shortest path and Danupon [START_REF] Patt-Shamir | Fast Routing Table Construction Using Small Messages[END_REF] give an algorithm with a smaller time complexity Õ(n (1/2) D (1/4) +D)-time, (1+o(1))approximation to solve the weighted single-pair shortest path computation. However, these algorithms are designed for a special topology or getting approximation results, which are less useful for logically centralized SDN routing environment.

3) SDN graph modeling: In SDN environments, a network can be represented as a graph. SDN controllers thus can compute the network information by utilizing the graph algorithms [START_REF] Raghavendra | Dynamic graph query primitives for sdn-based cloud network management[END_REF]. Most of the current SDN controllers such as POX, OpenDaylight and ONOS directly use Dijkstra or Floyd algorithm for path computation. Apparently, such designs lack scalability when processing all-pairs shortest path due to the O(n 3 ) time complexity.

4) Consistency model: The consistency problem is a notable trade-off in the distributed system. As for the distributed SDN controller, the consistency model can either be coarsegrained or fine-grained. Under coarse-grained design, consistency is guaranteed by controllers and is transparent to applications. For example, HyperFlow [START_REF] Tootoonchian | Hyperflow: a distributed control plane for openflow[END_REF] manipulates a publish/subscribe system and Kandoo [START_REF] Hassas Yeganeh | Kandoo: a framework for efficient and scalable offloading of control applications[END_REF] utilizes a hierarchical design (root controller and local controller) for maintaining the consistency of network information. Although ONOS [START_REF] Berde | ONOS: Towards an Open, Distributed SDN OS[END_REF] employs a Gossip protocol, which relaxes the constraints in data consistency of topology-related information, it cannot provide effective optimization for some particular operations like path computation. Moreover, ONOS delegates the complexity of maintaining the consistency to an external system while using an external system also increases the management difficulty. A fine-grained consistency policy is firmly correlated to a particular algorithm or application. A proper finegrained consistency model can increase the synchronization efficiency significantly. That is the reason why Onix [START_REF] Koponen | Onix: A Distributed Control Platform for Large-scale Production Networks[END_REF] does not synchronize network information in a transparent way for the SDN programmers, who have to create and maintain the network information by themselves explicitly.

VII. CONCLUSION

In this paper, we proposed ParaCon as a solution to address the performance bottleneck in the SDN control plane. To accelerate the path computation, we designed a set of mechanisms that utilize the untapped parallel computing potential of a multi-controller control-plane system. More specifically, ParaCon integrates the following features into an overall design: first, it exploits an asynchronous distributed algorithm for path computation in SDN; second, it reduces the transmission overhead by utilizing a hybrid consistent model. The evaluation results show that our design can achieve higher performance than legacy systems. Furthermore, we analyzed the convergence time, both theoretically and experimentally, comparing ParaCon with the state-of-the-art solutions. The results demonstrated the effectiveness of our architecture.

For future work, we are interested in integrating the design into commodity Open-Source SDN controllers, such as Open-Daylight or ONOS. Also, we plan to leverage the power of heterogeneous computing to further improve the performance.

  Definition 2: (Eventual Consistency Policy) After a control plane global information update H, supposed to change the status of the receiving controller from N O to Y ES, with an eventual consistency policy, the control plane is available also during the inconsistency window going from the instant when ∃!c ∈ C | p c = Y ES, and the instant when ∀c ∈ C, p c = Y ES.

Fig. 3 .

 3 Fig. 3. ParaCon structure. All controllers are distributed and sharing the global topology view. Switches are allocated to different controllers by an allocator, and each controller only computes the path information of its switches. All the controllers can compute the path together by parallel computing.

Fig. 4 .

 4 Fig. 4. The centralized version. There are several steps to compute the path.The right matrix shows the weight of shortest routes between switches. We use the blue square with the blue background to indicate that the path view is updated by the controller. After link 0, 2 has changed to 30, we put 0 and 2 to the queue Q. The algorithm gets the head of Q and computes the shortest route of it. After processing 0, we add the adjacent: 1 and 3 to Q because there are some changes when processing 0. If the queue is empty, the algorithm stops. In this example, the algorithm uses 6 steps: 0, 2, 1, 3, 0, 2 to finish the algorithm.

Algorithm 1 : 3 while 5 if now is a SYNC then 6 for 8 Send

 13568 Path computation (parallel, online, asynchronous) Input: Link view L, The number of switches n, The queue of switches need to be maintained Q id Output: Path view T 1 SYNC is a synchronization signal that makes the current controller synchronize the path view to other controllers. 2 while True do Queue Q id is not empty do 4 now ← Q id .head() Each controller i except current one id do 7 Send queue Q i to controller i. Path view T i, * to controller i.

9 continue 10 for

 910 Each adjacent node i of node now do 11 if T now,i can be relaxed then 12 Relaxing Path View T now,i . 13 Get t, which minimized prior calculation T now,i . 14 Forwarding Table F now,i ← t 15 if T now has changed then 16 Put all adjacent nodes k into Q Con(k) .

17

 17

Put

  SYNC into Q Con(k) .

1 Q0Fig. 5 .

 15 Fig. 5. The distributed version. 4 switches are associated with 2 controllers. We use blue color (blue queue, nodes, and matrix) to indicate controller C1 controls these elements; we also use the red color (red queue, nodes and matrix) to indicate controller C2 controls these elements. Thus, the matrix on the right is separated into two parts. The top two rows are maintained by controller C1 and the bottom two rows are maintained by controller C2. We use the square with the blue (red) background to indicate the path view is updated by controller C1 (C2). We also use the square with the purple background to indicate the path view is updated because of data synchronization from other controllers. For example, in step 2, '6' in C1 is synchronized from C2, and '35' in C2 is synchronized from C1. They only use 3 steps to finish the algorithm. Each controller computes the path and synchronizes the path information with other controllers at the end of each step.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Average computing time from 1 controller to 40 controllers. With the increasing of the number of switches, more controllers can get more benefit.

Fig. 9 .

 9 Fig.9. The correctness ratio in different topologies with 30 changes in the topology. The correctness ratio, which is defined as a percentage of consistent results, show the consistency influence on the result. The dot Changes in topology with green line indicate the time when our algorithm start processing the topology changes request from the queues in controllers, while controllers receive topology changes from switches and put them into queues in 1s interval. The initial process, before the first topology change, makes the correctness ratio start from 0% to 100%.

Fig. 10 .

 10 Fig. 10. The computing time per request with ONOS and ParaCon in different topologies and with 4 controllers

Jin

  Zhao received his B. Eng. Degree in computer communications from Nanjing University of Posts and Telecommunications, China, in 2001, and the Ph.D. Degree in computer science from Nanjing University, China, in 2006. He joined Fudan University in 2006. He stayed at University of Massachusetts Amherst for 1 year as a visiting scholar in 2014. His research interests include software defined networking, media streaming and network coding theory. He is a member of IEEE and ACM. Xin Wang received his Bachelor Degree and the Master Degree from Xidian University, Xian, China, in 1994 and 1997 respectively, in Information Theory and Communications. He received the Ph.D. Degree from Shizuoka University, Japan in 2002, in Computer Science. Since 2002, he has been with the School of Computer Science at Fudan University, where he is currently a full professor. Stefano Secci received the M.Sc. Degree in communications engineering from Politecnico di Milano, Milan, Italy, in 2005, and a dual Ph.D. Degree in computer science and networks from Politecnico di Milano and Telecom ParisTech, France. Now he is an Associate Professor at the Universite Pierre et Marie Curie (UPMC -Paris VI -Sorbonne Universites), Paris, France. Webpage: https://lip6.fr/ Stefano.Secci.

TABLE III THE

 III AMOUNT OF SYNCHRONIZED DATA IN 'CERNET' (WITH 30 TOPOLOGY CHANGES IN 1s INTERVAL, THE UNIT IS 1 MESSAGE DEFINED INSECTION II) Moreover, the amount of synchronized data in one controller is decreased with respect to the computation load, which is well balanced among all controllers. Furthermore, looking at the overhead due to one single change (the last row in the table), it is smaller than 1, 681 units (equal to n 2 with n = 41, which composes the whole path view) when the number of controllers equals to 4 and 6. The overhead is a bit larger than 1, 681 when the number of controllers equals to 2 (this overhead is introduced by Vector Clock and overwritten, as discussed in Section IV.C)

	Number of controllers	2	4	6
	Total data synchronized with 30 changes (global, all controllers)	129,601	149,117	193,069
	Average data synchronized with 30 changes (1 controller)	64,780	37,269	32,144
	Average data synchronized with 1 change (1 controller)	2,173	1,230	1,066
	in 1s interval, using the 'CERNET' topology. We use a number
	of controllers (2, 4 and 6) to perform our evaluation. Note
	that we balance switches (41 switches in 'CERNET') among
	all controllers. We export the number of messages that need
	to be synchronized in all controllers during path computation
	from our program to measure the amount of synchronized data.
	TABLE III reports the results, showing that the communication
	overhead due to synchronization marginally increases with an
	increasing number of controllers, knowing that one would not
	expect having a very high number of controllers in a common
	SDN scenario.

ACKNOWLEDGEMENTS

This work was partially funded by the 863 program (Grant no. 2015AA016106) and the EU FP7 IRSES MobileCloud Project (Grant no. 612212)


Gunther H Peichl
email: gunther.peichl@kfunigraz.ac.at

Rachid Touzani
email: rachid.touzani@univ-bpclermont.fr

AN ACCURATE FINITE ELEMENT METHOD FOR ELLIPTIC INTERFACE PROBLEMS

A finite element method for elliptic problems with discontinuous coefficients is presented. The discontinuity is assumed to take place along a closed smooth curve. The proposed method allows to deal with meshes that are not adapted to the discontinuity line. The (nonconforming) finite element space is enriched with local basis functions. We prove an optimal convergence rate in the H 1 -norm. Numerical tests confirm the theoretical results.

1. Introduction. Boundary value problems with discontinuous coefficients constitute a prototype of various problems in heat transfer and continuum mechanics where heterogeneous media are involved. The numerical solution of such problems requires much care since their solution does not generally enjoy enough smoothness properties required to obtain optimal convergence rates. Although fitted or adapted meshes can handle such difficulties, these solution strategies become expensive if the discontinuity front evolves with time or within an iterative process. Such a (weak) singularity appears also in the numerical solution of other types of problems like free boundary problems when they are formulated for a fixed mesh or for fictitious domain methods.

We address, in this paper, a new finite element approximation of a model elliptic transmission problem that allows nonfitted meshes. It is well known that the standard finite element approximation of such a problem does not converge with a first order rate in the H 1 -norm in the general case. We propose a method that converges optimally provided the interface curve is a sufficiently smooth curve. Our method is based on a local enrichment of the finite element space in the elements intersected by the interface. The local feature is ensured by the use of a hybrid approximation. A Lagrange multiplier enables to recover the conformity of the approximation. The derived method appears then rather as a local modification of the equations of interface elements than a modification of the linear system of equations. This property ensures that the structure of the matrix of the linear system is not affected by the enrichment.

Let us mention other authors who addressed this topic in the finite element context. We point out the so-called XFEM (eXtended Finite Element Methods) developed in Belytschko et al. [START_REF] Belytschko | Arbitrary discontinuities in finite elements, Internat[END_REF] where the finite element space is modified in interface elements by using the level set function associated to the interface. Such methods, that are used also for crack propagation, have in our point of view, the drawback of resulting in a variable matrix structure. Moreover, although no theoretical analysis is available, numerical experiments show that they are not optimal in terms of accuracy. Other authors like Hansbo et al. [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Hansbo | A Lagrange multiplier method for the finite element solution of elliptic interface problems using nonmatching meshes[END_REF], have similar approaches to ours but here also the proposed method seems to modify the matrix structure by enriching the finite element. In Lamichhane-Wohlmuth [START_REF] Lamichhane | Mortar finite elements for interface problems[END_REF] and Braess-Dahmen [START_REF] Braess | The Mortar element method revisited -What are the right norms ? Thirteen International Conference on Domain Decomposition Methods[END_REF], a similar Lagrange multiplier approach is used for a mortar finite element formulation of a domain decomposition method. Finally, in a work by Li et al [START_REF] Li | New cartesian grid methods for interface problems using the finite element formulation[END_REF], an immersed interface technique, inspired from finite difference schemes, is adapted to the finite el-ement context. Note also that the references where Lagrange multipliers are employed have for these multipliers as supports the edges defining the interface. In our method, the interface supports the added degrees of freedom but the Lagrange multipliers are defined on the edges intersected by the interface and thus serve to compensate the nonconformity of the finite element space rather than enforcing interface conditions, which are being naturally ensured by the variational formulation.

In the following, we use the space L 2 (Ω) equipped with the norm • 0,Ω and the Sobolev spaces H m (Ω) and W m,p (Ω) endowed with the norms • m,Ω and • m,p,Ω respectively. We shall also use the semi-norm

| • | 1,Ω of H 1 (Ω). Moreover, if Ω 1 and Ω 2 form a partition of Ω, i.e., Ω = Ω 1 ∪Ω 2 , Ω 1 ∩Ω 2 = ∅ and if v is a function in W m-1,p (Ω) with v |Ωi ∈ W m,p (Ω i ), then we shall adopt the convention v ∈ W m,p (Ω 1 ∪ Ω 2 ) and denote by v m,p,Ω1∪Ω2 the broken Sobolev norm v m,p,Ω1∪Ω2 = v m-1,p,Ω + v m,p,Ω1 + v m,p,Ω2 .
Similarly, we denote by • m,Ω1∪Ω2 and | • | m,Ω1∪Ω2 , the broken Sobolev norm and semi-norm respectively for the H m -space. Finally, we shall denote by C, C 1 , C 2 , . . . various generic constants that do not depend on mesh parameters and by |A| the Lebesgue measure of a set A and by A • the interior of a set A.

Let Ω denote a domain in R 2 with smooth boundary Γ and let γ stand for a closed

C 2 -curve in Ω which separates Ω into two disjoint subdomains Ω + , Ω -such that Ω = Ω + ∪ γ ∪ Ω -and ∂Ω + = γ. For given f ∈ L 2 (Ω) and a ∈ L ∞ (Ω) we consider the transmission problem:        -∇ • (a∇u) = f in Ω + ∪ Ω -, u = 0 on Γ, [u] = a ∂u ∂n = 0 on γ,
where [v] denotes the jump of a quantity v across the interface γ and n is the normal unit vector to γ pointing into Ω -. For definiteness we let

[v] = v --v + with v ± = v |Ω ± .
In addition to boundedness of the diffusion coefficient we assume

a ± ∈ W 1,∞ (Ω ± ), a(x) ≥ α > 0, for x ∈ Ω, (1.1)
i.e. a is uniformly continuous on Ω \ γ, but discontinuous across γ.

The standard variational formulation of this problem consists in seeking u ∈

H 1 0 (Ω) such that Ω a ∇u • ∇v dx = Ω f v dx ∀ v ∈ H 1 0 (Ω). (1.2)
In view of the ellipticity condition (1.1), Problem (1.2) has a unique solution u in H 1 0 (Ω) but clearly u / ∈ H 2 (Ω). We shall assume throughout this paper the regularity properties:

u |Ω -∈ H 2 (Ω -), u |Ω + ∈ H 2 (Ω + ), u 2,Ω -∪Ω + ≤ C f 0,Ω . (1.3)
Note that these assumptions are satisfied in the case where a |Ω -and a |Ω + are constants (see [START_REF] Lemrabet | Régularité de la solution d'un problème de transmission[END_REF][START_REF] Petzoldt | Regularity and error estimators for elliptic problems with discontinuous coefficients[END_REF] for instance).

In the following, we describe a fitted finite element method. defined by adding extra unknowns on the interface γ. It turns out that this method leads to an optimal convergence rate. Although it is well suited for the model problem it seems to be inefficient in more elaborate problems which, for example, involve moving interfaces. To circumvent this difficulty, we define a new method where the added degrees of freedom have local supports and then yield a nonconforming finite element method. We show that the use of a Lagrange multiplier removes this nonconformity and ensures an optimal convergence rate.

2. A fitted finite element method. Assume that the domain Ω is a convex polygon and consider a regular triangulation T h of Ω with closed triangles whose edges have lengths ≤ h. We assume that h is small enough so that for each triangle T ∈ T h only the following cases have to be considered:

1) T ∩ γ = ∅.

2) T ∩ γ is an edge or a vertex of T .

3) γ intersects two different edges of T in two distinct points different from the vertices. 4) γ intersects one edge and its opposite vertex. Let V h denote the lowest degree finite element space

V h = { v ∈ C 0 (Ω); v |T ∈ P 1 (T ) ∀ T ∈ T h , v = 0 on Γ},
where P 1 (T ) is the space of affine functions on T . A finite element approximation of (1.2) consists in computing

u h ∈ V h such that Ω a ∇u h • ∇v dx = Ω f v dx ∀ v ∈ V h . (2.1)
It is well known that, since u / ∈ H 2 (Ω), the classical error estimates (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF]) do not hold any more even though we still have the convergence result, lim h→0 uu h 1,Ω = 0.

A fitted treatment of the interface γ can however improve this result. Let for this purpose T γ h denote the set of triangles that intersect the interface γ corresponding to cases 3) and 4) above,

T γ h := {T ∈ T h ; γ ∩ T • = ∅},
and consider a continuous piecewise linear interpolation of γ, denoted by γ h , as shown in Figure 2.1. Clearly, γ h is the line that intersects γ at two edges of any triangle that contains γ. Unless the intersection of γ with the boundary of a triangle T does not coincide with an edge, T is split into two sets T + and T -separated by the curve γ. In case 3), the straight line γ h ∩ T splits T into a triangle K 1 and a quadrilateral that we split into two subtriangles K 2 and K 3 , where we choose K 2 such that K 1 ∩ K 2 = γ h . In case 4), γ h ∩ T splits T into two triangles K 1 and K 2 . In this case we set K 3 = ∅. This construction defines the new fitted finite element mesh of the domain Ω (see Figure 2.1). The splitting T = K 1 ∪ K 2 ∪ K 3 is not unique but the convergence analysis does not depend on it. Let us denote by T γ T the set of the three subtriangles of T . Below E h will stand for the set of all edges of elements and E γ h is the set of all edges that are intersected by γ (or γ h ), i.e.

E γ

h := {e ∈ E h ; γ ∩ e • = ∅}. For each T ∈ T h , E T is the set of the three edges of T . The fitted mesh is denoted by T F h , i.e.

K 1 K 3 K 2 γ h γ
T F h := T h ∪ T ∈T γ h ∪ K∈T γ K K .
and by S γ h := {T ; T ∈ T γ h }. Let us finally note that the curve γ h defines a new splitting of Ω into two subdomains Ω - h and Ω + h where Ω ± h is defined analogously to Ω ± with γ replaced by γ h .

Next we construct an approximation of the function a on the elements of T F h : For this purpose, let ã± be extensions of a ± to Ω such that ã± ∈ W 1,∞ (Ω). Such extensions exist due to the regularity of γ (see [START_REF] Adams | Sobolev Spaces[END_REF]). Define ãh ∈ W 1,∞ (Ω) by

ãh = ã+ in Ω + h , ã-in Ω - h ,
and denote by a h the piecewise linear interpolant of ã on T F h . Hence a h is continuous on Ω + h ∪ Ω - h and coincides with a on the nodes of T F h . In addition, the function a h is discontinuous across the line γ h and satisfies the properties,

a h|Ω + h ∈ W 1,∞ (Ω + h ), a h|Ω - h ∈ W 1,∞ (Ω - h ), (2.2)
a h 0,∞,Ω ≤ C a 0,∞,Ω , (2.3)
a h ≥ α > 0 a.e. in Ω. (2.4)
We now define the finite element space

W h = V h + X h , X h := {v ∈ C 0 (Ω); v |Ω\S γ h = 0, v |K ∈ P 1 (K) ∀ K ∈ T γ T , ∀ T ∈ T γ h }. Note that we have W h ⊂ H 1 0 (Ω).
A fitted finite element approximation is defined as the follows:

   Find u F h ∈ W h such that Ω a h ∇u F h • ∇v dx = Ω f v dx ∀ v ∈ W h . (2.5)
In order to study the convergence of Problem (2.5), we consider the auxiliary problem:

   Find u h ∈ H 1 0 (Ω) such that Ω a h ∇ u h • ∇v dx = Ω f v dx ∀ v ∈ H 1 0 (Ω).
(2.6)

We note that both problems (2.5) as well as (2.6) have a unique solution. The regularity properties (1.3) imply u + ∈ C 0 ( Ω+ ), u -∈ C 0 ( Ω-) and that u + and u - have a common trace on γ. Therefore u is continuous on Ω and the piecewise P 1 interpolant I h u ∈ W h is well defined. In the following let ũ± ∈ H 2 (Ω) stand for the extensions of u ± from Ω ± to Ω.

In the sequel, we assume that the fitted family of meshes (

T h ∪ T γ h ) h satisfies the condition h ̺ ≤ C h -θ (2.7)
for some θ ∈ [0, 1) and for which C is independent of h, where ̺ denotes the radius of the largest ball contained in any triangle in any triangle

T ∈ T F h . Lemma 2.1. Let u ∈ H 2 (Ω + ∪ Ω -).
1. We have the local interpolation error

|u -I h u| 1,T ≤ Ch |u| 2,T for T ∈ T h \ T γ h C h 2 ̺K (|ũ + | 2,K + |ũ -| 2,K ) for K ∈ T γ T , T ∈ T γ h , (2.8)
where ̺ K is the radius of the inscribed circle of K.

The global interpolation error is given by

|u -I h u| 1,Ω ≤ C h 1-θ |u| 2,Ω + ∪Ω -.
(2.9)

Moreover, if u ∈ W 2,∞ (Ω + ∪ Ω -) then |u -I h u| 1,Ω ≤ C h |u| 2,∞,Ω + ∪Ω -.
(2.10)

Proof. Since the local interpolation error estimate for T ∈ T h \ T γ h is classic in finite element theory (see [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF] or [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF] for instance), we only need to prove the second estimate on triangles where u is only piecewise smooth. Consider an element T ∈ T γ h and any subtriangle K ∈ T γ T . Without loss of generality we assume K ⊂ Ω + h , then

K = (K ∩ Ω + ) ∪ (K ∩ Ω -). Since K ∩ Ω -⊂ T ∩ Ω -∩ Ω + h and γ h interpolates the interface γ we obtain for the measure of K ∩ Ω - |K ∩ Ω -| ≤ |T ∩ Ω -∩ Ω + h | ≤ Ch 3 , (2.11)
with a constant C > 0 which depends on γ only. In view of I h u = I h ũ+ , the standard interpolation theory (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF] or [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]) implies

|u -I h u| 1,K ≤ |u -ũ+ | 1,K + |ũ + -I h ũ+ | 1,K ≤ |u -ũ+ | 1,K + C h 2 ̺ K |ũ + | 2,K .
(2.12)

Since ũ+ = u holds on K ∩ Ω + we obtain

|u -ũ+ | 1,K = |u -ũ+ | 1,K∩ Ω -≤ |u -| 1,K∩ Ω -+ |ũ + | 1,K∩ Ω -.
Applying Hölder's inequality with p = 3 2 and q = 3, the imbedding of H 1 (K) into L 6 (K) (Note that the imbedding constant can be bounded independently of h) and (2.11)

one can bound |u -| 1,K∩ Ω -(and analogously |ũ + | 1,K∩ Ω -) by |u -| 1,K∩ Ω -≤ |K ∩ Ω -| 1 3 ∇u - 0,6,K∩ Ω - ≤ C h ∇ũ - 0,6,K ≤ C h |ũ -| 2,K . Hence |u -ũ+ | 1,K ≤ C h (|ũ -| 2,K + |ũ + | 2,K ).
Inserting this estimate into (2.12) leads to

|u -I h u| 1,K ≤ C h 2 ̺ K (|ũ -| 2,K + |ũ + | 2,K ).
To prove the global interpolation error bound, we write

|u -I h u| 2 1,Ω = T ∈T h \T γ h |u -I h u| 2 1,T + T ∈T γ h K∈T γ T |u -I h u| 2 1,K ≤ Ch 2 T ∈T h \T γ h |u| 2 2,T + C T ∈T γ h K∈T γ T h 2 ̺ K (|ũ -| 2 2,K + |ũ + | 2 2,K ) ≤ C h 2 ̺ (|ũ -| 2 2,Ω + |ũ + | 2 2,Ω ) ≤ C h 2 ̺ |u| 2 2,Ω + ∪ Ω -,
where

̺ = min{̺ K : K ∈ T γ T , T ∈ T γ h }.

The calculation above indicates how the convergence rate can be improved in case

u ∈ W 2,∞ (Ω + ∪ Ω -) observing that |S γ h | ≤ Ch holds. Remark 2.1.
It is classic in finite element theory to assume that the meshes are regular in the sense that Condition (2.7) is satisfied for θ = 0. For the fitted meshes T γ h one cannot guarantee that such a condition is satisfied. To relax this constraint, we assume here (2.7) for a θ ∈ [0, 1) thus allowing a larger class of fitted meshes than permitted by θ = 0.

The following result gives the convergence rate for Problem (2.5). Theorem 2.1. Assume that the family of fitted meshes (T F h ) h satisfies the regularity property (2.7). Then we have the error estimate

|u -u F h | 1,Ω ≤ Ch 1-θ u 2,Ω + ∪ Ω -if u ∈ H 2 (Ω + ∪ Ω -), Ch u 2,∞,Ω + ∪ Ω - if u ∈ W 2,∞ (Ω + ∪ Ω -). (2.13)
Proof. We have from the triangle inequality

|u -u F h | 1,Ω ≤ |u -u h | 1,Ω + | u h -u F h | 1,Ω . (2.14)
To bound the first term on the right-hand side of (2.14), we proceed as follows: Let us subtract (2.6) from (1.2) and choose v = uu h . We have

Ω (a ∇u -a h ∇ u h ) • ∇(u -u h ) dx = 0.
Then

Ω a h |∇(u -u h )| 2 dx = - Ω (a -a h ) ∇u • ∇(u -u h ) dx = - Ω\S γ h (a -a h ) ∇u • ∇(u -u h ) dx - T ∈T γ h T (a -a h ) ∇u • ∇(u -u h ) dx.
The usual estimate for the interpolation error gives

a -a h 0,∞,Ω ≤ Ch ( ã 1,∞,Ω + h + ã 1,∞,Ω - h ) ≤ Ch a 1,∞,Ω + ∪Ω -.
with a constant C which only depends on a reference triangle, (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF], p. 124). Thus we obtain

Ω\S γ h (a-a h ) ∇u•∇(u-u h ) dx ≤ C h a 1,∞,Ω + ∪Ω -|u| 1,Ω\S γ h |u-u h | 1,Ω\S γ h . (2.15)
Next we consider a triangle T ∈ T γ h which we split as

T = (T ∩ Ω + ∩ Ω + h ) ∪ (T ∩ Ω -∩ Ω - h ) ∪ (T ∩ Ω + ∩ Ω - h ) ∪ (T ∩ Ω -∩ Ω + h ).
As before, we obtain

T ∩Ω + ∩Ω + h (a-a h ) ∇u•∇(u-u h ) dx ≤ Ch a 1,∞,Ω + ∪Ω -|u| 1,T ∩Ω + ∩Ω + h |u-u h | 1,T ∩Ω + ∩Ω + h .
Arguing as in the proof of Lemma 2.1, the generalized Hölder inequality together with (2.11) yields the estimate

T ∩Ω + ∩Ω - h (a -a h ) ∇u • ∇(u -u h ) dx = T ∩Ω + ∩Ω - h (a + -a - h )∇u + • ∇(u + -u h ) dx ≤ C a 0,∞,Ω |T ∩ Ω + ∩ Ω - h | 1/3 ∇u + 0,6,T ∩Ω + ∩Ω - h ∇(u + -u h ) 0,T ∩Ω + ∩Ω - h ≤ C h a 0,∞,Ω |ũ + | 2,T ∇(u + -u h ) 0,T .
Analogous estimates hold with + andinterchanged. Collecting the four contributions to the triangle T one obtains

T (a -a h ) ∇u • ∇(u -u h ) dx ≤ Ch ( a 0,∞,Ω + a 1,∞,Ω + ∪ Ω -) × |ũ + | 2,T ∇(u + -u h ) 0,T ∩Ω + + |ũ -| 2,T ∇(u --u h ) 0,T ∩Ω -.
Combining this estimate with (2.15) leads to

Ω a h |∇(u -u h )| 2 dx ≤ Ch a 1,∞,Ω + ∪ Ω -|u| 1,Ω\S γ h |u -u h | 1,Ω\S γ h + C h ( a 0,∞,Ω + a 1,∞,Ω + ∪ Ω -) × T ∈T γ h |ũ + | 2,T ∇(u + -u h ) 0,T ∩Ω + + |ũ -| 2,T ∇(u --u h ) 0,T ∩Ω - ≤ C h a 1,∞,Ω + ∪ Ω -|u| 1,Ω\S γ h |u -u h | 1,Ω\S γ h + C h ( a 0,∞,Ω + a 1,∞,Ω + ∪ Ω -)(|ũ + | 2,S γ h + |ũ -| 2,S γ h ) ∇(u -u h ) 0,S γ h ≤ C h ( a 0,∞,Ω + a 1,∞,Ω + ∪ Ω -) |u| 2,Ω + ∪ Ω -∇(u -u h ) 0,Ω , which by (2.4) implies |u -u h | 1,Ω ≤ C h ( a 0,∞,Ω + a 1,∞,Ω + ∪ Ω -) |u| 2,Ω + ∪ Ω -. (2.16) To bound the norm | u h -u F h | 1
,Ω , we have from problems (2.6) and (2.5),

Ω a h ∇( u h -u F h ) • ∇v dx = 0 ∀ v ∈ W h .
Standard finite element approximation theory combined with (2.2)-(2.3) gives

| u h -u F h | 1,Ω ≤ C inf v∈W h | u h -v| 1,Ω , (2.17)
which together with (2.16) implies

| u h -u F h | 1,Ω ≤ C | u h -I h u| 1,Ω ≤ C | u h -u| 1,Ω + C |u -I h u| 1,Ω ≤ C h |u| 2,Ω + ∪ Ω -+ C |u -I h u| 1,Ω .
The interpolation error is bounded using (2.9) or (2.10).

A hybrid approximation.

The method presented in the previous section has proven its efficiency as numerical tests will show in the last section. In more elaborate problems like time dependent or nonlinear problems where the interface γ is a moving front, the subtriangulation T γ h moves within iterations and then the matrix structure has to be frequently modified. To remedy to this difficulty, we resort to a hybridization of the added unknowns. More specifically, the added discrete space X h is replaced by a nonconforming approximation space. In addition, a Lagrange multiplier is used to compensate this inconsistency. The hybridization enables to locally eliminate the added unknowns in each triangle T ∈ T γ h . In the sequel we fix an orientation for the interface γ. This induces an orientation of the normals to the edges e ∈ E γ h by following the interface in the positive direction. The jump of a function v across an edge e ∈ E γ h can then be defined as

[v] e (x) := lim s→0,s>0 v(x + sn(x)) -lim s→0,s<0 v(x + sn(x)) ≡ v + (x) -v -(x), x ∈ e,
where n is the unit normal to e.

To develop this method, we start by defining an ad-hoc formulation for the solution u h of (2.6). Let us define the spaces

Z h := H 1 0 (Ω) + Y h , Y h := {v ∈ L 2 (Ω); v |Ω\S γ h = 0, v |T ∈ H 1 (T ) ∀ T ∈ T γ h , [v] = 0 on e, ∀ e ∈ E h \ E γ h }, Q h := e∈E γ h H -1 2 00 (e),
where

H -1 2
00 (e) is the dual space of the trace space

H 1 2 00 (e) := {v |e ; v ∈ H 1 (T ), e ∈ E T , v = 0 on d ∀ d ∈ E T , d = e}.
We remark that the jumps [v] for v ∈ Z h can be interpreted in H

1 2
00 (e) for e ∈ E γ h . This is due to the fact that v ∈ H 1 (T ) for all T ∈ T h , that for every e ∈ E γ h , the jump of v lies in H 1 2 (e) and vanishes at the endpoints of e as well as on at least two adjacent edges. This motivates the choice of Q h .

The elements of Q h will be referred to by µ = (µ e ) e∈E γ h . We endow Z h with the broken norm

u Z h = ( T ∈T h |u| 2 1,T ) 1/2 .
On Q h we use the norm

µ Q h = e∈E γ h µ e 2 H -1 2
00 (e) .

Above, the integrals over edges e are to be interpreted as duality pairings between H 00 (e). We mention that the broken norm in Z h reflects the fact that Z h is not a subspace of H 1 0 (Ω). Next we define the variational problem,

Find ( u H h , λ h ) ∈ Z h × Q h such that: T ∈T h T a h ∇ u H h • ∇v dx - e∈E γ h e λ h [v] ds = Ω f v dx ∀ v ∈ Z h , (3.1)
e∈E γ h e µ [ u H h ] ds = 0 ∀ µ ∈ Q h . (3.2)
The saddle point problem (

u H h Z h + λ h Q h ≤ C f 0,Ω , (3.3)
with a constant C which is independent of h. Proof. Problem (3.1)-(3.2) can be put in the standard variational form

A ( u H h , v) + B(v, λ h ) = (f, v) ∀ v ∈ Z h , B( u H h , µ) = 0 ∀ µ ∈ Q h , where
A (u, v) = T ∈T h T a h ∇u • ∇v dx, B(v, µ) = - e∈E γ h e µ [v] ds, (f, v) = Ω f v dx.
The bilinear form A is clearly continuous and coercive on the space Z h × Z h . The bilinear form B is also continuous on

Z h × Q h .
Next we verify that B satisfies the inf-sup condition, i.e. there exists δ > 0 such that for every λ ∈ Q h there exists v µ ∈ Z h such that Let χ T denote the characteristic function of T and define

B(v µ , µ) ≥ δ v µ Z h µ Q h i.e. e∈E γ h e µ e [v µ ] ds ≥ δ v µ Z h µ Q h (3.
v µ = T ∈T Γ h χ T v T .
Since there are as many edges in E γ h as triangles in T γ h then [v µ ] = v T holds for every edge e ∈ E γ h . Hence we obtain

µ 2 Q h = e∈E γ h µ e 2 -1/2,e ≤ T ∈T γ h ∇v T 2 0,T = e∈E γ h e µ e [v µ ] ds.
Furthermore,

v µ 2 Z h = T ∈T h ∇v T 2 0,T = T ∈T γ h ∇v T 2 0,T
holds. This implies

µ 2 Q h v µ 2 Z h ≤ T ∈T γ h ∇v T 2 0,T 2 = B(v µ , µ) 2 .
Adjusting the sign of v µ this is equivalent to (3.4) We are now able to present a numerical method to solve the interface problem. This one is simply derived as a finite element method to solve the saddle point problem (3.1)-(3.2). We consider for this end a piecewise constant approximation of the Lagrange multiplier. Let us define the finite dimensional spaces,

Z h := V h + Y h , Y h := {v ∈ L 2 (Ω); v |Ω\S γ h = 0, v |K ∈ P 1 (K) ∀ K ∈ T γ T , ∀ T ∈ T γ h , [v] = 0 on e, ∀ e ∈ E h \ E γ h }, Q h := µ ∈ e∈E γ h L 2 (e); µ |e = const. ∀ e ∈ E γ h .
The hybrid finite element approximation is given by the following problem:

Find (u H h , λ h ) ∈ Z h × Q h such that: T ∈T h T a h ∇u H h • ∇v dx - e∈E γ h e λ h [v] ds = Ω f v dx ∀ v ∈ Z h , (3.6
)

e∈E γ h e µ [u H h ] ds = 0 ∀ µ ∈ Q h . (3.7)
Let us give some additional remarks before proving convergence properties of this method.

1. The matrix formulation of the method has the following form

  A C 0 C T D B 0 B T 0     u v λ   =   b c 0   , (3.8)
where the vector u contains the values of u H h at nodes of the mesh T h , i.e. components of u H h in the Lagrange basis of V h , v contains the components of u F h in the basis of Y h , and λ has as components the values of λ h on the edges of E γ h . There is clearly no simple method to eliminate off diagonal blocks in the system (3.8) in order to decouple the variables. More specifically, our aim is to eliminate the unknowns v.

2. The method must be viewed in the context of an iterative process like the Uzawa method, where the Lagrange multiplier λ h is decoupled from the primal variable u H h . In such situations, each iteration step consists in solving an elliptic problem with a given λ h . Let us recall that, due to the local feature of the basis functions of nodes on edges of E γ h , the unknowns associated to these nodes can be eliminated at the element level. This is a basic issue in our method.

3. We point out that equation (3.7) entails

[u H h ] = 0 on e, ∀ e ∈ T γ h . (3.9)
This follows from the fact that u H h is an affine function on each edge of T γ h . This implies that actually

u H h ∈ W h . Choosing v ∈ W h in (3.1) we find Ω a h ∇u H h • ∇v dx = Ω f v dx.
This yields u H h = u F h . 4. Convergence analysis. This section is devoted to the proof of existence, uniqueness and stability of the solution of (3.6)-(3.7) as well as its convergence to Problem (3.1)- (3.2).

For this result we need a localized quasi-uniformity of the mesh. More precisely, we assume that

|e| ≥ Ch ∀ e ∈ E γ h . (4.1)
In addition, we make the following assumption:

The distance of the intersection point of γ with any edge e ∈ E γ h to the endpoints of e can be bounded from below by δh, where δ is independent of h.

(4.2)
Although this assumption appears to be quite restrictive, numerical tests have shown that it can be actually ignored in applications. Theorem 4.1. Assume that the family of meshes (T h ) h satisfies Property (4.1). Then Problem (3.6)-(3.7) has a unique solution. Moreover, we have the bound

u H h Z h + λ h Q h ≤ C f 0,Ω , (4.3)
where the constant C is independent of h. Proof. It is clearly sufficient to prove the inf-sup condition (see for instance Brezzi-Fortin [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]):

sup v h ∈Z h \{0} e∈E γ h e µ h [v h ] ds v h Z h µ h Q h ≥ β > 0 ∀ µ h ∈ Q h . (4.4)
In the following, for each triangle T ∈ T γ h , we shall denote by e + T (resp. e - T ) the edge where γ enters T (resp. leaves T ), and by ẽT the remaining edge of T (see Figure (4.1)). Recall that we fixed an orientation for γ. The trace inequality (see [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF], eq. (2.5)) and the Poincaré inequality owing to v = 0 on ẽT , yield for

T ∈ T γ h , v 0,e ± T ≤ C 4 h -1 2 v 0,T + h 1 2 ∇v 0,T ≤ C 5 h 1 2 ∇v 0,T . (4.7)
On the other hand, Assumption (4.2) implies the uniform boundedness of ∇ϕ e ± T 0,T . From (4.6) and (4.7) we obtain then ∇v h 0,T ≤ C 6 ∇v 0,T .

Using the inf-sup condition (3.4) and (4.5), we finally obtain

µ h Q h v h Z h ≤ C 6 µ h Q h v Z h ≤ C 7 e∈E γ h e µ h [v] ds = C 7 e∈E γ h e µ h [v h ] ds.
Finally, obtaining the estimate (4.3) is a classical task that we skip here.

We now prove the main convergence result. Theorem 4.2. Assume hypotheses (2.7) and (4.1) are satisfied, then there exists a constant C, independent of h, such that

u -u H h Z h ≤ Ch 1-θ |u| 2,Ω + ∪Ω -if u ∈ H 2 (Ω + ∪ Ω -), Ch u 2,∞,Ω + ∪Ω -if u ∈ W 2,∞ (Ω + ∪ Ω -).
Proof. From classical theory of saddle point problems (see [START_REF] Girault | Finite Element Methods for the Navier-Stokes Equations, Theory and Algorithms[END_REF], p. 114), we obtain from Theorem 4.1,

u H h -u H h Z h + λ h -λ h Q h ≤ C inf v∈Z h u H h -v Z h + inf µ∈Q h λ h -µ Q h . (4.8)
Furthermore, using (Braess [4], Theorem 4.8), Property (3.9) implies that Estimate (4.8) can be improved, for the error on u H h by

u H h -u H h Z h ≤ C inf v∈Z h u H h -v Z h . (4.9)
To bound the right-hand side, we choose v = I h u, where I h is the previously defined Lagrange interpolant in Z h . Since u H h = u h (see Theorem 3.1), then by using (2.9) and (2.16), This can be achieved without assuming (4.2). In this case, no error estimate is to be expected for the Lagrange multiplier.

u H h -I h u Z h = u h -I h u Z h ≤ u -I h u Z h + u -u h Z h ≤ C 1 h 1-θ |u| 2,Ω + ∪Ω -+ C 2 h |u| 2,Ω + ∪Ω -. If u ∈ W 2,∞ (Ω + ∪ Ω -), then Estimate (2.10) yields u H h -I h u Z h ≤ C h u 2,∞,Ω + ∪Ω -. ( 4
Finally, since the Lagrange multiplier λ h can be interpreted in terms of u h (see Theorem 3.1), it is interesting to see how good is its approximation λ h . Let, for this, E h denote the set

E h := e∈E γ h e.
Theorem 4.3. Under the same hypotheses as in Theorem 4.2, we have the following error bounds

λ h -λ h Q h ≤    C(h 1-θ + h 1 2 ) |u| 2,Ω + ∪Ω -+ λ 0,E h if λ h ∈ L 2 (E h ), Ch 1-θ |u| 2,Ω + ∪Ω -+ λ 1 2 ,E h if λ h ∈ H 1 2 (E h ).
Proof. We use the abstract error bound (4.8). Let, for e ∈ E γ h , λ e := 1 |e| e λ h ds.

Using Lemma 7 in Girault-Glowinski [START_REF] Girault | Error Analysis of a Fictitious Domain Method Applied to a Dirichlet Problem[END_REF], we obtain the bound

λ h -λ e H - 1 2 00 (e)
≤ Ch

1 2 λ h 0,e if λ h ∈ L 2 (e), and
λ h -λ e H -1 2
00 (e)

≤ Ch λ h 1 2 ,e if λ h ∈ H 1 2 (e),
Combining these bounds with (4.8), (4.9) and (4.10) achieves the proof.

A numerical test.

To test the efficiency and accuracy of our method, we present in this section a numerical test. We consider an exact radial solution and test convergence rates in various norms.

Let Ω denote the square Ω = (-1, 1) 2 and let the function a be given by a

(x) = α if |x| < R 1 , β if |x| ≥ R 1 ,
where α, β > 0. We test the exact solution

u(x) =      1 4α (R 2 1 -|x| 2 ) + 1 4β (R 2 2 -R 2 1 ) if |x| < R 1 , 1 4β (R 2 2 -|x| 2 ) if |x| ≥ R 1 .
We choose R 1 = 0.5 and R 2 = √ 2. The function f and Dirichlet boundary conditions are determined according to this choice. Note that unlike the presented model problem, we deal here with non homogeneous boundary conditions but this cannot affect the obtained results.

The finite element mesh is made of 2N 2 equal triangles. According to the definition of a, the interface γ is given by the circle of center 0 and radius R 1 . The error is measured in the following discrete norms:

e 0,h := 1 M M i=1 (u(x i ) -u h (x i )) 2 1 2 , e 0,∞ := max 1≤i≤M |u(x i ) -u h (x i )|, e 1,h := T ∈T h T |I h (∇u)(x) -∇u h | 2 1 2
, where x i are the mesh nodes, M is the total number of nodes, and I h is the piecewise linear interpolant. We denote in the sequel by p the ratio α/β. Table 1 presents convergence rates for the standard P 1 finite element method using the unfitted mesh (2.1) with the choice p = 1/10.

h 2.02 × 10 -2 6.28 × 10 -2 20 6.78 × 10 -3 1.05 1.09 × 10 -2 0.89 5.23 × 10 -2 0.26 40 3.61 × 10 -3 0.91 5.81 × 10 -3 0.91 3.68 × 10 -2 0.51 80 1.83 × 10 -3 0.98 3.06 × 10 -3 0.92 2.56 × 10 -2 0.52 160 9.44 × 10 -4 0.95 1.55 × 10 -3 0.98 1.82 × 10 -2 0.49 Table 1. Convergence rates for a standard (unfitted) finite element method.

As expected, numerical experiments show poor convergence behavior. Let us consider now the results obtained by the present method, i.e. (3.1)-(3.2) or equivalently (2.5). We obtain for p = 1/10 and p = 1/100 the convergence rates illustrated in Tables 1 and2 2.2 5.22 × 10 -5 1.7 3.50 × 10 -4 1.5

Table 3. Convergence rates for the hybrid finite element method with p = 1/100.

Tables 2 and3 show convergence rates that are even better than the theoretical results. This is probably due to the choice of a discrete norm but may also be due to a superconvergence phenomenon. Rates for the L 2 -norm give also good behavior. For the L ∞ -convergence rate, we can note that we dot retrieve the second order obtained for a continuous coefficient problem. However, these rates are better (1.5 rather than 1) than the ones obtained for a standard finite element method and, moreover, the error values are significantly lower in our case. It is in addition remarkable that the error values depend very weekly on p but the convergence rates are independent of this value.

6. Concluding remarks. We have presented an optimal rate finite element method to solve interface problems with unfitted meshes. The main advantage of the method is that the added unknowns that deal with the interface singularity do not modify the matrix structure. This feature enables using the method in more complex situations like in problems with moving interfaces. The price to pay for this is the use of a Lagrange multiplier that adds an unknown on each edge that cuts the interface. This drawback can be easily removed by using an iterative method such as the classical Uzawa method or more elaborate methods like the Conjugate Gradient. The good properties of the obtained saddle point problem enable choosing among a wide variety of dedicated methods. This topic will be addressed in a future work. Let us also mention that the present finite element method does not specifically address problems with large jumps in the coefficients. These ones are in addition ill conditioned and this drawback is not removed by this technique.

Fig. 2 . 1 .

 21 Fig. 2.1. Subdivision of interface triangles.

|∇v T | 2 dx, which implies µ e 2 - 1 / 2

 212 µ = (µ e ) e∈E γ h ∈ Q h and an edge e ∈ E γ h choose a triangle T ∈ T γ h which has e as one of its edges. Define v T ∈ H 1 (T ) as the solution of  dx = e µ e ϕ ds for ϕ ∈ H 1 e (T ) where H 1 e (T ) = {ϕ ∈ H 1 (T ); ϕ = 0 on ∂T \ e}. By Green's theorem we obtain µ e -1/2,e = ∂v T ∂n -1/2,e ≤ ∇v T 0,T , e µ e v T ds = T ,e ≤ T |∇v T | 2 dx = e µ e v T ds.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

σFig. 4 . 1 .( 4 . 5 )

 4145 Fig. 4.1. Definition of e + T , e - T , and ẽT .Let µ h ∈ Q h , and let v ∈ Z h be the function given by Problem (3.5). We define a function v h ∈ Z h by   

  3.1)-(3.2) indicates that the continuity of u h across the edges of E γ h is enforced by a Lagrange multiplier technique. Theorem 3.1. Problem (3.1)-(3.2) has a unique solution ( u H h , λ h ) ∈ Z h × Q h . Moreover, we have u H h = u h and the following estimate holds

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . e + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . T γ h . T -T + e + γ h e -e +
	e -	e -

  .10)

	Remark 4.1. As it was previously mentioned, we know that if Problem (3.1)-
	(3.2) has a unique solution (u H h , λ h ) then u H h = u F h where u F h is the solution of Problem (2.5) and therefore the error estimate (2.13) holds. Consequently, Theorem 4.2 can be
	simply proven by obtaining a nonuniform inf-sup condition (i.e. (4.4) with β = β(h)).

Table 2 .

 2 respectively. Convergence rates for the hybrid finite element method with p = 1/10.

	h -1	e 0,h	Rate	e 0,∞	Rate	e 1,h	Rate
	10 3.45 × 10 -3 20 8.18 × 10 -4 40 1.70 × 10 -4 80 3.94 × 10 -5 160 8.57 × 10 -6	2.1 2.3 2.1 2.2	4.25 × 10 -3 1.72 × 10 -3 5.22 × 10 -4 1.64 × 10 -4 4.89 × 10 -5	1.3 1.7 1.7 1.7	1.75 × 10 -2 6.87 × 10 -3 2.81 × 10 -3 1.02 × 10 -3 3.59 × 10 -4	1.3 1.3 1.5 1.5
	h -1	e 0,h	Rate	e 0,∞	Rate	e 1,h	Rate
	10 3.26 × 10 -3 20 7.91 × 10 -4 40 1.72 × 10 -4 80 4.01 × 10 -5 160 8.82 × 10 -6	2.0 2.2 2.1	4.07 × 10 -3 1.74 × 10 -3 5.47 × 10 -4 1.74 × 10 -4	1.2 1.7 1.6	1.69 × 10 -2 6.65 × 10 -3 2.72 × 10 -3 9.88 × 10 -4	1.3 1.3 1.5

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .