Keywords: Multi-Hop Wireless Networks, Low-Power and Lossy Networks, Energy Efficiency, TCP

Low-power and Lossy Networks (LLNs), like wireless networks based upon the IEEE 802.15.4 standard, have strong energy constraints, and are moreover subject to frequent transmission errors, not only due to congestion but also to collisions and to radio channel conditions.

This paper introduces an analytical model to compute the total energy consumption in an LLN due to the TCP protocol. The model allows us to highlight some tradeoffs as regards the choice of the TCP maximum segment size, of the Forward Error Correction (FEC) redundancy ratio, and of the number of link-layer retransmissions, in order to minimize the total energy consumption.

List of Figures

Energy consumption with short or long TCP segments, as a function of the number of link layer attempts r (with B = 5 × 10 -4 ). . . . . . . . . . 6

Consumed energy using short or long TCP segment, as a function of the redundancy ratio α (B = 3 × 10 -4 , h = 5). . . . . . . . . . . . . . . . . . 7

Long (MSS=512 bytes) versus short (MSS=64 bytes) in a multi-hop TCP transmission: prefer the short MSS above the curves, the long one below. 8

Long (MSS=512 bytes) versus short (MSS=64 bytes) in a multi-hop TCP transmission: prefer the short MSS above the curves, the long one below.

List of Tables

Introduction

IP-based Low-power and Lossy Networks (LLNs) are composed of (possibly many) nodes which are subject to strong constraints on power and computing capabilities. Besides, link-layer technologies like IEEE 802.15.4 [START_REF][END_REF], often used for building wireless LLNs, impose very low limits on link frame sizes. Many usual applications for LLNs (e.g., wireless sensor networks) do not always require full reliability. However, there exist several use-cases for LLNs where the reliability of data delivery can be critical, like over-the-air software updates [START_REF] Dunkels | Distributed TCP Caching for Wireless Sensor Networks[END_REF][START_REF] Park | GARUDA: Achieving Effective Reliability for Downstream Communication in Wireless Sensor Networks[END_REF][START_REF] Wan | Pump-slowly, fetch-quickly (PSFQ): a reliable transport protocol for sensor networks[END_REF]. In such cases, the Transmission Control Protocol (TCP) [START_REF] Postel | Transmission Control Protocol[END_REF] is a common choice for ensuring an end-to-end reliable transport in IP-based LLNs.

TCP is the most used transport protocol in IP networks. It provides a reliable transfer service for all kinds of applications. TCP can manage data that is damaged, duplicated, lost, or delivered out of order. This is achieved by assigning a sequence number to each octet transmitted and requiring a positive acknowledgement (ACK) from the recipient. If the TCP ACK segment is not received within the timeout interval, the data segment is retransmitted.

This paper focuses on the energy cost of reliability when TCP is used in multi-hop wireless LLNs. In what follows, we will illustrate some main issues by considering the case of TCP in IPv6-enabled LLNs based on the 6LoWPAN protocols [START_REF] Kushalnagar | IPv6 over Low-Power Wireless Personal Area Networks (6LoWPANs): Overview, Assumptions, Problem Statement, and Goals[END_REF].

Many parameters affect the energy consumption of a reliable transport protocol, like channel conditions (e.g., wireless losses, collisions), the level of link-layer reliability, the number of links/hops, or the maximum size of link-layer frames (e.g., 127 bytes for IEEE 802.15.4). The length of such frames may require the fragmentation of IPv6 datagrams before sending them. The IETF 6LoWPAN working group has defined a new protocol layer between the IPv6 and the medium access control (MAC) layers. The 6LoWPAN layer compresses the IPv6 header and fragments IPv6 packets into short MAC frames. So, if the size of a TCP segment exceeds the maximum allowed length, the 6LoWPAN layer fragments it to make it fit into link-layer frames. Therefore, sending short TCP segments does not require any fragmentation. However, the use of a small maximum segment size1 (MSS) increases the number of TCP segments and, thus, the number of TCP acknowledgements and the corresponding MAC frames sent. Besides, with small MSS values the overhead due to TCP headers becomes larger; this, coupled with the small size of MAC frames, may result in a very low protocol efficiency.

This paper introduces a mathematical model aimed at predicting the energy consumed by the wireless nodes of an LLN in a bulk-data transfer scenario, with TCP used as a reliable transport layer. The model estimates TCP energy performance based on the bit error rate, the maximum number of retransmissions at the link layer, the number of hops between the sender and the receiver, the amount of FEC, and the TCP maximum segment size.

The proposed model allows us to study the tradeoffs involved in sending short versus long TCP segments. We assume that the energy consumed in a data transfer depends mainly on the number of bits sent. Thus, our analytical model estimates the number of bits sent by all nodes in the network, taking into account the cumulated cost of all link-layer transmissions. Indeed, the number of MAC frames sent can be large with respect to the initial amount of data to transmit, because of link-layer retransmissions, the redundancy added for Forward Error Correction (FEC), and also due to end-to-end (i.e., transport layer) retransmissions when a TCP segment is lost before reaching the receiving node.

We apply the model to study the energy efficiency of TCP over an LLN using 6LoW-PAN and 802.15.4 protocols, and study the effect of the TCP segment size, of the FEC redundancy ratio, and of the maximum link layer retransmission on the total energy consumption.

The remainder of the paper is organized as follows. Section 2 gives a brief overview of the related work in the area of TCP performance in LLNs, and in multihop wireless networks in general. Section 3 presents the derivation of our analytical model. In Section 4, we apply the analytical results to derive the best TCP segment size strategy in terms of energy consumption. Finally, in Section 5, we give our conclusions and provide directions for future work.

Related work

There exist many studies of the performance of TCP in multi-hop wireless networks (see for example [START_REF] Hanbali | A Survey of TCP over Ad Hoc Networks[END_REF] and references therein), though many of them are concerned mostly with parameters like TCP throughput, and with protocol enhancements aimed at improving such performance metrics. Another amount of work focuses on the computational cost of TCP (e.g., [START_REF] Wang | Computational Energy Cost of TCP[END_REF]), or on the comparison of the energy efficiency for different TCP congestion control algorithms (like SACK, Tahoe, Reno and NewReno) over wireless networks (e.g., [START_REF] Heng | Energy Efficiency Comparison of SACK, Tahoe, Reno and NewReno Over Ad Hoc Network Based on Analytic Hierarchy Process[END_REF][START_REF] Seddik-Ghaleb | A performance study of TCP variants in terms of energy consumption and average goodput within a static ad hoc environment[END_REF][START_REF] Zorzi | Energy efficiency of TCP in a local wireless environment[END_REF]). In [START_REF] Fu | The Impact of Multihop Wireless Channel on TCP Performance[END_REF], Fu et al. study in the impact of the MAC layer protocol on the TCP throughput over wireless multi-hop networks. The authors show that, in contrast to wired networks, the TCP throughput does not increase by increasing the window size, and that (under the conditions of the study) better performance is achieved with a window size equal to the number of hops divided by four.

In [START_REF] Lilakiatsakun | Enhancing TCP energy efficiency for mobile hosts[END_REF], Lilakiatsakum and Senevirane propose an energy-efficiency metric to compare the performance of different versions of TCP. The metric is the ratio between the amount of bits sent by all nodes and the size of the application data. In this work, we adopt a variant of that metric (the total amount of bits sent); we use a fixed value for the application data size in all our numerical and simulation scenarios.

Bansal et al. [START_REF] Bansal | Energy efficiency and capacity for TCP traffic in multi-hop wireless networks[END_REF] propose an analytical model to compute the energy consumed for carrying a TCP flow over a multi-hop wireless network. The authors assume that all wireless links have the same packet error rate P and the same transmission power. They compute the energy spent by all nodes for sending a single TCP segment. The energy consumption is then a function of the packet error rate, the number of hops, and the maximum number of link-layer retransmissions. However, [START_REF] Bansal | Energy efficiency and capacity for TCP traffic in multi-hop wireless networks[END_REF] does not take into account the cost of sending link-layer acknowledgements, nor the cost of transport-layer acknowledgements, which we add in the model introduced here. Besides, we explicitly consider the impact of lower-layer fragmentation and of error correction, as well as different per-link error rates on the energy-efficiency of TCP.

Barman et al. [START_REF] Barman | TCP Optimization through FEC, ARQ, and Transmission Power Tradeoffs[END_REF] and Gallucio et al. [START_REF] Galluccio | An Analytical Study of a Tradeoff Between Transmission Power and FEC for TCP Optimization in Wireless Networks[END_REF] present an analytical study of a TCP optimiza-tion problem in an hybrid wired/wireless network where the last hop is a wireless link. The authors of both papers define an utility function which is the ratio of the throughput to the cost of a TCP connection. Our work completes these studies by introducing a multi-hop model for computing the energy consumption. Note that in this paper we are not interested in TCP throughput, because data rates are a secondary concern in LLNs; instead, we focus mainly on the energy costs of a TCP connection over multiple lossy links.

TCP energy consumption model

This section introduces a mathematical model to estimate the energy consumed by a TCP transmission in a wireless LLN. We assume that such energy mainly corresponds to data emission and reception, and thus directly depends on the number of bits sent by all nodes.

We therefore compute the expected total amount of bits sent for a successful end-toend TCP data transmission, as a function of several network parameters, namely: the bit error rate, the link-layer maximum number of attempts, the FEC redundancy ratio, the number of hops between the source and receiver TCP hosts, and the TCP maximum segment size.

For the convenience of the reader, Table 1 lists most of the variables used in this paper.

Link layer: one-hop model

We first focus on modeling the link-layer energy consumption, considering a CSMA-CA network with error correction control techniques combining Automatic Repeat Request (ARQ) [START_REF] Fairhurst | Advice to link designers on link Automatic Repeat reQuest (ARQ)[END_REF] and Forward Error Correction (FEC).

Link layer mechanisms

Standard ARQ uses the cyclic redundancy check (CRC) error-detecting code that is added to the data: The receiver uses the error-detecting code number to check the integrity of the received data. After receiving a correct frame, the receiver replies by an ACK. If the sender does not receive any ACK before a timeout2 -because either the original message or the ACK is lost, or they contain errors-, the sender retransmits again the same message. If the receiver sees that the frame is damaged the receiver discards it and does not send an ACK. The ARQ algorithm continues until the sender receives an ACK or exceeds a predefined number of attempts r.

ARQ is the algorithm most frequently used by link-layer protocols to reduce the packet error rate. However, if the wireless network becomes very lossy, ARQ would increase the transmission delay between the source and the receiver. Abrupt increases of the end-toend delay increase the round-trip time and may lead to a spurious TCP timeout. This can deteriorate the TCP performance. Expected number of bits sent for a successful end-to-end packet transmission knowing that it has succeeded S Average number of bits sent for successfully transmitting a TCP segment An orthogonal approach consists in applying Forward Error Correction (FEC). FEC is a good solution for decreasing the packet error rate. The main idea of FEC is to add redundancy to the original frame, to allow the destination node to detect and correct some bit errors. In our case, if the size of a network datagram is greater than the maximum transmission unit (MTU) of the link layer, the datagram is divided into fragments of length K bits, and the FEC algorithm adds (α × K) redundancy bits to form a frame of length D. The ratio α = D-K K between the amount of redundancy added by FEC and the original frame length is called the redundancy ratio.

In what follows we will consider an error-correction method like the well-known Reed-Solomon [START_REF] Reed | Polynomial codes over certain finite fields[END_REF] (RS) algorithm. By adding D -K bits to the K data, the RS algorithm can correct up to (D -K)/2 bits.

Performance of one-hop transmissions

First, we consider here only a transmission between two immediate neighbors, without intermediate nodes. In CSMA-CA, two types of messages are used during data transmission: the data message (i.e., the message that contains the useful data), and the acknowledgement message sent by the receiving node. Figure 1 shows the three possible cases: a. Failure: A failure due to the loss of the data frame. The sender will retransmit the frame.

b. Partial failure: The data frame is correctly received by the receiver, while the acknowledgement frame is lost 3 . Therefore, the sender will (needlessly) retransmit the data frame.

c. Success: A successful transmission of both the data and the acknowledgement frames.

Sender Receiver Da ta fra m e . The probability of a successful transmission of a link-layer frame depends on the bit error rate, the size of the message, and the redundancy ratio. Here we assume that all bits undergo independent random errors with a given bit error rate B, each bit being correctly received with a probability 1 -B. Moreover, we assume that a frame transmission fails as soon as the number of bits erroneously received exceeds the number of correctable bits. The probability P fail of having a failure (case (a) in Fig. 1) is then the error probability for a data frame, i.e.,

P fail := 1 - c i=0 D i B i (1 -B) D-i
where D is the number of bits of a data frame, c = D-K 2 = αK 2 is the number of correctable bits, and for two integers i and j, j i is the number of possibilities of choosing i elements among j, i.e.,

j i = j! i!(j -i)! .
Likewise, denoting by A the size of an acknowledgement frame, the probability of a partial failure is the probability that the data frame is correctly received but the acknowledgement frame contains errors:

P partial := (1 -P fail )(1 -(1 -B) A ),
whereas the probability of a success is

P succ := (1 -P fail )(1 -B) A .
Note that we assume that acknowledgement frames are sent without adding redundancy. This is a reasonable assumption, given the current practices in most wireless technologies.

In order to estimate the energy consumption, we compute the expected number of bits sent at the link layer. Remark that in cases of success or partial failure, D + A bits are sent, whereas only D bits are sent in cases of failure. Recall however that the redundancy ratio determines the number of useful bits per frame, and thus the number of frames to send.

We now take into account the link-layer sending attempts, following the ARQ technique implemented by the MAC protocol. We denote by r the maximum number of attempts, after which the sender of a data frame considers the receiver is unreachable.

The reception failure probability of the data frame after r attempts is simply the probability F of having r successive failures, i.e.:

F := P r fail . (1)
In that case, the total number of bits sent only comes from the sender (no acknowledgement is sent), and thus equals

H f := r × D. (2)
With probability 1 -F , the receiver gets the data frame within the r link-layer attempts. In that case, the total number of bits sent depends on the number of failures and partial failures before a success, if any (we only know here that at least one attempt led to a success or a partial failure). There are two possibilities:

• either all r attempts were failures or partial failures, but at least one was a partial failure (since we are in the case where the receiver got the data),

• or the last of the k ≤ r attempts was a success (in the sense of case (c) in Fig. 1).

Conditioned on the receiver getting the data frame within the r link-layer attempts, the expected total number of bits sent H s can therefore be computed as follows:

H s = 1 1 -F ( r i=1 r i P i partial P r-i fail (rD+iA)+ r k=1 P succ k-1 i=0 k -1 i P i partial P k-1-i fail (kD+(i+1)A)).
(3) In (3), k stands for the index of the first success in r attempts (if any), and i is the number of partial failures among all attempts.

Multi-hop model

Let us now focus on the multi-hop case. An end-to-end transmission succeeds if the message reaches the destination after a certain number h of hops.

In this paper, we assume that link layer transmissions on each hop are independent. We denote B i the bit error rate and F i the frame error rate of the i th hop. F i is computed as per (1), taking B = B i . Therefore, the probability that a frame is correctly received by a destination node is simply the probability Q s that all h one-hop transmissions succeed, i.e.,

Q s = h i=1 (1 -F i ) ( 4
)
where h is the number of hops from the sender to the receiver. We assume here that the MAC layer is able to detect duplicate frames, such as those that are produced in case of a partial failure (case (b) in Fig. 1). This can be implemented, for instance, by using a sequence number in the MAC frame headers; when a node receives from one neighbor two successive frames with the same sequence number, it assumes that the corresponding Ack frame was lost and deletes the second frame 4 . This avoids the propagation of several copies of the same data frame.

Figure 2 shows the two possibilities for the outcome of the end-to-end transmission of a frame. Sender Receiver . . .

(a) End-to-end failure scenario: the frame cannot be forwarded after r unsuccessful retransmissions.

Sender

Receiver .

(b) End-to-end success scenario: the frame arrives at the destination. This scenario may also include partial failures over one or more hops (not depicted). Again, we express the expected number of bits sent, conditioned on the success of the end-to-end transmission (encompassing possible link-layer retransmissions, within the limit of r total attempts per hop).

• Knowing that the destination node correctly receives the message, the expected number of bits sent by all network nodes is simply

E s := h i=1 H s i . (5)
• Knowing that the message was lost in one of the h hops, the number of bits sent depends on the hop where the loss (i.e., the failure of all r attempts) occurs. The expected value then equals

E f := h k=1 ( k-1 i=1 H s i + H f k ) k-1 j=1 (1 -F j )F k 1 -Q s ( 6
)
where H s i and H f i are computed using ( 3) and (2), respectively.

TCP performance

We now study the energy consumption of TCP in a multi-hop wireless network, taking into account the fact that TCP segments may be fragmented by lower layers if the total size of data frames (that depends on the TCP maximum segment size (MSS)) exceeds the MTU of those layers.

We can intuitively think of several opposite effects of the MSS, which suggest that a trade-off has to be found: i) If TCP segments are fragmented, the loss of one fragment leads to the loss of the original TCP segment and therefore a new TCP end-to-end retransmission. Further, using short TCP segments (i.e., a small MSS) saves CPU power associated to fragmentation and reassembly. Moreover, sending several fragments (treated independently by lower layers) increases the probability of collision between two TCP data fragments, due to the hiddennode problem. As an illustration of that phenomenon, one can refer to Fig. 4, provided later on, where the number of collisions with long TCP segments is much larger than when short TCP segments are used. ii) On the other hand, the use of long TCP segments reduces the number of TCP segments sent by the source, and thus, the TCP header overhead (which can be very large, for small MTU values) and the number of TCP acknowledgements.

To analyze such trade-off, we now apply the model presented in Section 3.2. Our objective here is to determine an MSS choice that optimizes TCP performance, in terms of energy consumption.

Let us consider that the MSS of TCP and the MTU of the MAC layer are such that each TCP segment is fragmented into m link-layer frames. This encompasses the particular case when no fragmentation occurs, that simply corresponds to m = 1. We will make the assumption that TCP ACKs always fit in a single MAC frame, though it would be easy to consider a more general case in which TCP ACKs are also fragmented. Note finally that we do not take into account the use of the delayed ACKs mechanism by the TCP receiver; including the effect of such a mechanism is left to future work.

To simplify the analysis, we assume that all fragments of a TCP data segment have the same total size D DATA , and all TCP acknowledgement frames have the same total size D ACK . At the destination, the segment is reconstructed from the received fragments. The loss of a single link-layer frame induces the loss of the whole TCP segment, and thus the retransmission of all m fragments. Further, we consider that the number of TCP retransmissions is not limited; that is, the TCP source keeps on sending a segment until it receives a TCP acknowledgement from the TCP receiver.

From the previous analysis, the success probability P s of a TCP segment transmission attempt is simply the probability that all m data fragments be correctly sent to the destination, and the TCP ACK be successfully sent back to the source:

P s = Q m s × Q s,ack
, where Q s and Q s,ack denote the success probability of the multi-hop transmission of a TCP data fragment and of a TCP acknowledgement frame, respectively. Q s and Q s,ack are simply obtained by applying ( 4), but replacing D by respectively D DATA and D ACK .

In our model, each TCP data fragment transmission succeeds or fails independently of the others. Knowing that a transmission is successful at the TCP level (i.e., the TCP ACK is correctly received by the TCP source, which implies that all m fragments correctly reached the destination), the expected total number of bits sent by all nodes equals:

S s := E s × m + E s,ack
As above, E s and E s,ack are obtained from (5) using D DATA and D ACK , respectively, as the size of a frame. Likewise, we also define E f,ack from [START_REF] Fu | The Impact of Multihop Wireless Channel on TCP Performance[END_REF].

Knowing that a TCP transmission attempt has failed, the expected number of bits sent end-to-end by all nodes is:

S f := 1 1 -P s I f (1 -Q m s )
end-to-end transmission failure of one or more of the m fragments

+ (E s × m + E f,ack )Q m s (1 -Q s,ack )
end-to-end transmission failure of the TCP ACK

, where I f is the expected total number of bits sent for the m fragments to reach the destination, knowing that they (i.e., at least one) finally fail. Formally,

I f := m k=1 m k kE f + (m -k)E s (1 -Q s ) k (Q s ) m-k .
After some algebra, we get

I f = m(1 -Q s )E f + mE s Q s (1 -Q m s ).
To simplify the analysis, we will assume that the TCP window is equal to one TCP segment. This assumption is justified because a small window is a sensible choice for networks with a moderate number of hops [START_REF] Fu | The Impact of Multihop Wireless Channel on TCP Performance[END_REF]. Moreover, such small windows are typical of current TCP implementations for LLNs (e.g., Contiki OS [START_REF] Dunkels | Contiki -a lightweight and flexible operating system for tiny networked sensors[END_REF]), since the memory and CPU constraints of wireless embedded devices make it difficult to fully implement TCP's congestion control mechanisms.

We can now compute the total number of bits that have to be sent by all nodes, to successfully transmit both a TCP segment and its corresponding TCP ACK. Since we assumed the number of TCP retransmissions is unbounded (i.e., the TCP sender keeps on retrying until the transmission succeeds), the mean number of transmissions for a given TCP segment equals 1/P s . This therefore corresponds to a total number of bits sent (per segment) of S := S f (1/P s -1) + S s .

Finally, to successfully send a given amount M of application ("useful") data, the expected number of bits sent by all wireless nodes is the number of segments M/MSS , multiplied by the expected number of bits sent for a successful segment transmission. It is that final value that will be considered as representative of the overall energy consumption of the TCP transmission.

Results and discussion

We will now compare the predictions of our analytical model to simulation results, and discuss the tradeoff between sending long or short segments in different scenarios.

We apply our model to study the energy consumption of TCP on IPv6 over Low-Power Wireless Personal Area Networks [START_REF] Kushalnagar | IPv6 over Low-Power Wireless Personal Area Networks (6LoWPANs): Overview, Assumptions, Problem Statement, and Goals[END_REF], considering IEEE 802.15.4 [START_REF][END_REF] as the link-layer technology. The maximum link-layer frame size is thus equal to 127 bytes. The 6LoWPAN layer adapts the IPv6 datagrams to the link-layer MTU, as explained in [START_REF] Kushalnagar | IPv6 over Low-Power Wireless Personal Area Networks (6LoWPANs): Overview, Assumptions, Problem Statement, and Goals[END_REF].

We consider two MSS choices (MSS=64 and MSS=512 bytes) for a given TCP session. When the MSS is 64 bytes, no fragmentation is performed by the 6LoWPAN layer, whereas in the case where MSS=512 bytes, this adaptation layer splits each segment into 8 frames.

Energy model

To compute the energy consumed by wireless nodes to successfully send a data frame, we specify in this section a simple energy model. We denote by n the average number of neighbors of each node5 . If a frame is sent by a wireless node, it is received by its n neighbors (assumed thus to be in an active listening state). Therefore the energy consumed to send a link-layer frame equals

E := (TransmitEnergy + n × ReceiveEnergy) × FrameSize
Table 2 shows the energy consumption values that we considered in our numerical computations. Those values correspond to a transmit and receive consumption currents of

Model assessment

We first validate the results of our analytical model, through the simulation of the following scenario. A TCP sender sends a short file (51.2 Kbytes) to a TCP receiver. We consider an average size of 20 bytes for IPv6 headers, thanks to the LOWPAN IPHC [START_REF] Montenegro | Transmission of IPv6 Packets over IEEE[END_REF] compression. The simulation results plotted are average values after 30 simulation runs.

The TCP window size is set to 1, which implies that no congestion losses occur (we only consider one flow here). As a result, retransmissions can only be due to transmission errors, as described in Section 3. For the sake of simplicity, all links have the same bit-error rate B i = B, i = 1, . . . , h. Unless indicated otherwise, the parameters used in simulations and in numerical computations correspond to the default values shown in Table 3. Finally, as in Section 3.3, the maximum number of transport-layer retransmissions is not bounded. Figure 3 shows that analytical results closely match simulations results when a short TCP segment size is chosen. However, we can see that with long TCP segments, the model tends to underestimate the energy consumption. This difference comes from the presence of collisions between fragments of a given TCP data segment, that are not taken into account in the analytical model but do occur in simulation, as illustrated in Fig. 4. In the figure, nodes 1 and 6 are the source and destination TCP nodes, respectively, and node 2 to 5 are intermediate nodes.

Moreover, Fig. 3 shows that using short TCP segments becomes interesting, from an energy point of view, when the bit-error rate is high (above 10 -4 ). For example, for a BER of 4 × 10 -4 , the total consumed energy with MSS = 512 bytes is three times larger than with MSS = 64 bytes. For low values of B, the number of retransmissions is small, hence the energy consumption becomes roughly independent of B. We remark on Figs. 3 and5 that simulation results fit better our analytical results when no fragmentation is performed by the link layer. This is again due to collisions that occur in the large-MSS case and are not encompassed by our model. Fig. 5 also illustrates that increasing the number of link-layer attempts decreases the energy consumption, especially so when there is fragmentation (i.e., a large MSS). Indeed, giving the link layer more chances to send a data frame to its next-hop node reduces the discard probability of TCP segments, and therefore the number of end-to-end retransmissions.

In the following, we study the impact of the other system parameters on the energy consumption, based on the model only. We should therefore stay aware that the energy consumption may be a little underestimated for large values of the MSS.

FEC redundancy ratio and energy consumption

We first study the impact of FEC on the energy consumption. Fig. 6 shows that there seems to be an optimal amount of redundancy, in terms of energy consumption. Below the optimal value, adding redundancy reduces the probability of losses, and thus reduces the energy consumption. When α is above such optimal value, however, energy expenditure steadily increases because of the redundancy overhead.

Recall that the MTU being fixed, the number of frames to send depends on the redundancy ratio according to a stairstep function, hence the discontinuities in the figure.

Selecting the TCP MSS to minimize energy consumption

Depending of the numerous parameters of a given scenario, it appears that the same MSS value is not always the most efficient one in terms of energy. We now intend to summarize the effect of all parameters, by focusing on the best MSS strategy to implement. In the following figures, we compare the two values M SS = 64 and M SS = 512, and we concentrate on the boundary values, that delimitate zones where one MSS value outperforms the other. We plot those frontier curves in Figs. 7 and8. In each figure, the area above the curve represents the case where a TCP MSS value of 64 bytes consumes less energy than an MSS of 512 bytes, while the opposite holds below the curves. Figure 7 shows the two zones, depending on the transmission distance and the BER. We remark that for a given BER, short MSSs tend to outperform long MSSs when the distance grows: it is more and more interesting to use short MSS values instead of long ones. Indeed, for large networks the cost of end-to-end retransmissions will exceed the potential economy in terms of TCP header overhead. For a given distance, short MSSs are better suited for large BER environments since (energy-spending) segment retransmissions tend to occur more frequently. As previously observed, increasing the maximum number of link layer attempts r reduces the one-hop transmission failures, and thus limit the effects of errors, hence favoring long MSSs.

Finally, Fig. 8 illustrates the effect of FEC mechanisms. Not surprisingly (the FEC reducing the effect of transmission errors), redundancy makes large MSSs outperforms small MSSs due to the overhead reduction they allow.

Conclusions and Future Work

In this paper, we have proposed an analytical model to estimate the number of bits sent by all wireless nodes in a TCP session in a Low power and Lossy Network, in order to evaluate the overall energy consumption. The model has been validated through simulations, using the INETMANET framework [START_REF] Vagas | INETMANET framework[END_REF] of the OMNet++ network simulator, in the context of TCP over 6LoWPANs.

Our main outcomes regard the choice of an energy-saving Maximum Segment Size (MSS) for TCP. We have shown that using a large TCP segment size is less energyconsuming in small, low-error networks, while it becomes interesting to reduce the MSS when the network is large or very lossy. The impact of the number of attempts at the link layer, as well as the use of FEC, has also been studied.

The work presented here could be extended in several ways. A first interesting direction would be to model the collision process, that has been observed in our simulations for large MSS values. We would also like to consider the case when duplicate frames are not detected at the link layer; likewise, the transport layer modeling could be extended to encompass TCP's delayed acknowledgement mechanism, and also larger TCP windows. Finally, we are currently investigating an adaptation algorithm which dynamically adjusts the TCP segment size to the optimal MSS value.

  (a) Failure. Sender Receiver D a t a fr a m e Ack fra me . (b) Partial failure. Sender Receiver D a t a fr a m e A ck fr am e (c) Success.

Figure 1 :

 1 Figure 1: Failure and success scenarios for one link-layer transmission attempt.

Figure 2 :

 2 Figure 2: Failure and success scenarios in a multi-hop transmission.

Figure 3 :

 3 Figure 3: Energy consumption with long or short TCP segments, as a function of the BER B.

Figure 4 :

 4 Figure 4: Number of collisions in a multi-hop scenario.

Figure 5 :

 5 Figure 5: Energy consumption with short or long TCP segments, as a function of the number of link layer attempts r (with B = 5 × 10 -4 ).

Figure 6 :

 6 Figure 6: Consumed energy using short or long TCP segment, as a function of the redundancy ratio α (B = 3 × 10 -4 , h = 5).

Figure 7 :

 7 Figure 7: Long (MSS=512 bytes) versus short (MSS=64 bytes) in a multi-hop TCP transmission: prefer the short MSS above the curves, the long one below.

Figure 8 :

 8 Figure 8: Long (MSS=512 bytes) versus short (MSS=64 bytes) in a multi-hop TCP transmission: prefer the short MSS above the curves, the long one below.

1

  Failure and success scenarios for one link-layer transmission attempt. . . 2 Failure and success scenarios in a multi-hop transmission. . . . . . . . . . 3 Energy consumption with long or short TCP segments, as a function of the BER B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Number of collisions in a multi-hop scenario. . . . . . . . . . . . . . . . . 5

Table 1 :

 1 Notations used in this paper; capital italics letters correspond to probabilities, bold letters to (expected) numbers of bits.

	Variable Definition
	h	Number of hops between source and destination
	r	Maximum number of transmission attempts at the link layer
	m	Number of fragments corresponding to a single TCP segment (due
		to link layer fragmentation)
	α	FEC redundancy ratio
	B	Bit error rate
	P fail	Probability of a failure in a link-layer transmission attempt of a
		data frame
	P partial	Probability of a link-layer data frame being correctly received and
		the corresponding acknowledgement being lost
	P succ	Probability of a successful link-layer transmission attempt
		(data+acknowledgement are correctly received)

F Probability that a destination node does not receive a link layer data frame after r attempts Q s Probability of an end-to-end packet transmission success Q f Probability of an end-to-end packet transmission failure D Link-layer data frame size A Link-layer acknowledgement frame size H f Expected number of bits sent after r attempts knowing that the (one-hop) transmission has failed H s Expected number of bits sent within r attempts knowing that the (one-hop) transmission has succeeded E f Expected number of bits sent for an end-to-end packet transmission knowing that it has failed E s

Table 2 :

 2 Energy parameters.

	Parameter	Value
	Transmit Energy 0.24 µJ/bit
	Receive Energy	0.21 µJ/bit
	20 mA and 17.7 mA, respectively, a tension of 3V and a transmission rate of 250 kbit/s
	(thus, a per-bit transmission duration of 4 µs). Such figures are typical of low-power
	wireless modules currently on the market (e.g., Crossbow TelosB 6 ) .

Table 3 :

 3 Default simulation parameters

	Parameter	Value
	h	5
	r	3
	α	0
	BER B	3 × 10 -4
	Link-layer Ack frame size	40 bits
	Link-layer data frame header 120 bits
	IP header	160 bits
	TCP header	160 bits

The TCP MSS option refers to the maximum amount of data a TCP segment can carry, i.e., excluding the TCP header and any options that might be present.

Remark that we do not consider time issues in this model, thus, only losses can lead to retransmissions.

In a multi-hop setting, the receiver would nonetheless relay the frame to the next hop.

Extending the current model to consider the case where there is no duplicate frame detection is left as future work. Note that the current 802.15.4 specification does not explicitly mandate such a procedure, though it is possible to implement it.

In the plotted curves, we considered n = 2.

http://www.hoskin.qc.ca/uploadpdf/Instrumentation/CrossBow/hoskin TPR2400CA 42efb73715b8b.pdf

Acknowledgements

The work of A. Ayadi has been funded by the Pôle de Recherche Avancée en Communications (PRACom). We would like to thank Ahmed Triki for his help with the simulations.


MD Jean-Alexandre Long
email: jalong@chu-grenoble.fr

PhD Nikolai Hungr

PhD Michael Baumann

MD Jean-Luc Descotes

MD Michel Bolla

PhD Jean-Yves Giraud

MD Jean-Jacques Rambeaud

PhD Jocelyne Troccaz

à la diffusion de documents scientifiques

Keywords: Brachytherapy, Focal Therapy, Image fusion, Motion Tracking, Prostate cancer, Robotics, prostate biopsies

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

Development of a Novel Robot for Transperineal Needle Based Interventions: Focal Therapy, Brachytherapy and

Prostate Biopsies.

INTRODUCTION

Numerous methods exist for treating localized prostate cancer. Among these, radical prostatectomy, external beam radiotherapy and prostate brachytherapy are considered standard treatments, while focal therapies such as focal cryotherapy and photodynamic therapy, although more recent, have been shown to have the potential to improve functional outcomes after treatment of a localized cancer [START_REF] Eggener | Focal therapy for prostate cancer: possibilities and limitations[END_REF][START_REF] Lazzeri | Focal therapy meets prostate cancer[END_REF] . Except HIFU, most of these treatments use implanted needles. Brachytherapy and ablative treatments use a template to insert the needles along a grid of horizontal holes, the depth of each needle being adjusted visually using 2D trans-rectal ultrasound (TRUS) guidance. Limitations to this implantation technique may exist. Both the insertion of the needles and the movement of the probe cause significant motion and deformation of the prostate [START_REF] Stone | Prostate gland motion and deformation caused by needle placement during brachytherapy[END_REF] , which could cause implantation inaccuracies. Needle insertion is also restricted to the template's horizontal axis, thereby restricting access behind the pubic arch and limiting eligibility to patients with prostates < 60 cm3 [START_REF] Heidenreich | EAU guidelines on prostate cancer. Part 1: screening, diagnosis, and treatment of clinically localised disease[END_REF][START_REF] Thompson | Guideline for the management of clinically localized prostate cancer: 2007 update[END_REF] .

We propose a computer-assisted robotic system called PROSPER, that consists of a robotic needle insertion device and a static 3D ultrasound probe [START_REF] Hungr | Design of an ultrasound-guided robotic brachytherapy needle-insertion system[END_REF][START_REF] Reynier | MRI/TRUS data fusion for prostate brachytherapy. Preliminary results[END_REF][START_REF] Baumann | Towards 3D ultrasound image based soft tissue tracking: a transrectal ultrasound prostate image alignment system[END_REF] (Figure 1). The robot allows needles to be inserted whatever the implanted angles, prostate size and prostate mobility, hence improving implantation accuracy, allowing for targeting of tumor foci and extending the indications for larger prostates. Although focal therapy is the long term goal for the system, brachytherapy was chosen as the first application because of the current lack of long term data concerning focal therapy [START_REF] Arvold | Risk of death from prostate cancer after radical prostatectomy or brachytherapy in men with low or intermediate risk disease[END_REF] . This device could also be available for multiparametric prostate biopsies.

This paper presents the robot characteristics and procedure. We report the results of tests to assess the accuracy of needle insertion in prostate phantoms.

MATERIALS AND METHODS

ROBOT DESIGN

The robotic needle manipulator consists of a needle holder that is mounted on the lateral side of the US probe. It positions the needle along the appropriate insertion axis and drives it to a given depth. The robot uses seven motors, giving it seven degrees of freedom (figure 2). To increase safety, a mechanical release system disengages the needle from the robot in case the latter comes in contact with the pubic bone. The needle guide at the front of the needle insertion module and the bushing holding the needle are sterilizable.

The rest of the robot is covered by a sterile plastic cap. The robot has been designed and developed by our group and described in detail in a previous paper [START_REF] Hungr | Design of an ultrasound-guided robotic brachytherapy needle-insertion system[END_REF] .

CLINICAL PROCEDURE DESCRIPTION

An Ultrasonix™ 4DEC9-5 end-fire 3D probe is connected to an Ultrasonix RP™ machine.

The robotic needle manipulator is rigidly calibrated pre-operatively to the probe placed in the rectum. A 3D US reference volume is acquired. The needle trajectories and seed positions are then defined with respect to the reference prostate, according to the expected planning.

The robot positions the needle-holder in front of the perineum and inserts the first needle. In case of pubic arch interference, the needle is withdrawn and a partial re-planning is done to modify the needle trajectory in order to avoid the pubic arch. Once the needle has been inserted to its planned position, a verification procedure is applied to check for any prostate motion or deformation caused by the insertion. This is done by taking a second 3D US volume and registering it with the reference volume. If the target is not reached, the needle depth is adjusted iteratively until the expected point is reached. The treatment is applied, and the procedure is repeated for the rest of the needles.

PHANTOM TESTS

We performed a set of simulated procedures by depositing glass beads through inserted needles into realistic, mobile synthetic prostate phantoms. The results were measured by segmenting the inserted beads in CT scan volumes of the phantoms. Glass beads were chosen instead of titanium seeds because of the absence of artifacts in the CT images.

We developed the phantoms to obtain realism in US, CT and MRI images, but also in the mechanical soft-tissue behavior during needle insertion. The phantoms were built of soft PVC plastic, and included a hard perineum, soft peri-prostatic region, a harder prostate with an echogenic capsule and a hollow rectum for probe insertion. The differences in elasticity of the various anatomical regions made for realistic needle insertion forces and prostate motions. The prostate was molded with ten 1mm diameter glass beads embedded inside, acting as targets for needle insertion. The US probe was inserted into the phantom's rectum, giving realistic images of the prostate, with the glass targets being easily segmented. The phantom and its targets could also be seen clearly in CT images, allowing us to verify placement in CT volumes (figure 3). A detailed description of the phantoms has been published elsewhere [START_REF] Hungr | A realistic deformable prostate phantom for multi-modal imaging and needleinsertion procedures[END_REF] .

The experiment conducted with the PROSPER system involved trying to insert needles as close as possible to 90 targets embedded in 9 phantoms. As this study is an exploratory trial, no validated information on the effect of size and on the variability of needle accuracy was available, thus the sample size had no statistical rationale. We planned initially to construct 10 phantoms, each with 10 target beads embedded inside. Unfortunately, given the complexity of building an anthropomorphic synthetic prostate phantom, one phantom was damaged and was not available for the study.

The experimental procedure was as follows: the targets were segmented in a reference 3D US image and their coordinates were sent to the robot, which proceeded to insert the respective needle accordingly. An 18 gauge Mick Ripple-Hub™ needle was used, with an insertion speed of 5mm/s and a rotation speed of 8rps.

Needle rotation is known to decrease needle-tissue forces, needle bending and tissue deformation [START_REF] Meltsner | Observations on rotating needle insertions using a brachytherapy robot[END_REF] . Once the needle was inserted, a second US volume was acquired and the initial reference volume registered to it.

The new deformed target location was sent back to the robot. Then, it adjusted the needle depth to the closest point along the needle axis to this new target location. This was repeated until no further depth change could be made. A glass bead was then deposited through the hollow needle to mark the needle's final position.

Targets that were near the superior surface of the prostate were approached at a 10˚ horizontal and vertical

inclination

RESULTS

In the phantom experiments, at the apex, the median distance between the centers of the target and inserted beads was 2.28mm, while at the base it was 3.83mm (p<0.001). The more the targets were on the left, the higher the error was (2.2 vs 3.6 mm, p=0.002). The median amount by which the target moved in the needle insertion direction after image registration was 4.0 mm at the apex and 6.5 mm at the base of the prostate (median 5.46 mm, IQR 3.77-6.71). This can be understood as the amount of prostate motion in the needle-insertion direction. A large majority of the targets required one single depth correction to reach the final insertion point, while about 5% required two or more corrections. Accuracy was not different between horizontal and angled needle inclinations (2.7 vs 2.82 mm, p=0.28) and between anterior and posterior locations (2.7 vs 2.7 mm, p= 0.59). The maximum corrections provided by the tracking system were found in left (7.1 mm), base (6.5 mm) and posterior locations (5.9 mm) (table 1). Homogeneous motions of the prostate were detected by the tracking algorithm in all three spatial directions, x, y and z. This suggested, not just prostate translation in the needle direction (z), but also prostate rotation, that was accentuated with the depth of the needle (punctures at the base provide a bigger motion than those at the apex, whatever the axis considered) (table 2).

DISCUSSION

An important aspect for successful transperineal needle-based interventions is the ability to track prostate motion intra-operatively. Deformations and motions are highly underestimated in conventional practice.

Artignan et al. estimated that the prostate moved by at least 5 and 3 mm respectively in a sagittal and coronal plane, related only to the probe insertion [START_REF] Artignan | Online ultrasound image guidance for radiotherapy of prostate cancer: impact of image acquisition on prostate displacement[END_REF] . Moreover, according to Roy, needle insertion could result in motions of the base of the prostate of up to 15 mm [START_REF] Roy | Prostate gland motion, deformation and distorsion during needle placement in patients undergoing permanent seed implantation: Implications for pre-planning[END_REF] . It has been stated in the literature that the primary axis of prostate mobility during an insertion is along the needle insertion axis [START_REF] Taschereau | Seed misplacement and stabilizing needles in transperineal permanent prostate implants[END_REF] . Our experiments confirmed that prostate rotation can also affect the results: the deeper is the insertion, the more the prostate rotated. It was also noticed that peripheral needles caused more rotation than central needles, as could be expected, if the prostate is assumed to be constrained by the pubo-prostatic ligaments.

Clinical outcomes of this imprecision may be negligible in a treatment of the whole gland like for brachytherapy, but they may be more serious in focal therapy or targeted biopsies.

The results described above show that in our synthetic phantoms, needle insertion caused significant motion, on the order of 5mm in the direction of the needle axis. Without the registration step, the resulting seed distribution would have been significantly offset from the planned distribution. By correcting the needle depth based on the prostate motion, the accuracy of the system in the needle insertion direction was measured to be 2.28 mm at the apex and a little more (3.83 mm) at the base due to poorer image resolution at greater depths. The experiments showed that the system was capable of drastically reducing the errors caused by prostate mobility in the cranial-caudal direction. This was true for the different depths of insertion and approach angles tested.

We defined optimal needle insertion and rotation speed according to a needle insertion force study previously published [START_REF] Hungr | Design of an ultrasound-guided robotic brachytherapy needle-insertion system[END_REF] . In the phantom (similarly to in vivo force measurements reported in the literature during live brachytherapy) [START_REF] Podder | Vivo motion and force measurement of surgical needle intervention during prostate brachytherapy[END_REF] , the maximal force (1.8N) was reached when traversing the perineum. The amount of prostate motion was minimal with an insertion speed of 5mm/s. Although this speed is lower than during a manual procedure, it allows a safe insertion. Indeed, the insertion can be interrupted at any time by pressing the emergency stop button.

The reason why the robot is less precise when targeting to the left of the prostate is unclear. A possible explanation is that the weight of the robot might create increased gravitational strain on the robot causing a curvature towards the left side and thus a lower accuracy. This will have to be improved in the upcoming version of the robot.

The accuracy results are hardly comparable with systems previously described. Indeed, PROSPER is the only robot providing a motion tracking system and angulated needle insertion. As a result it could improve pubic collision avoidance, providing a better access to anterior and lateral parts of the prostate and extending the indications for larger prostates [START_REF] Song | Robotic needle guide for prostate brachytherapy: clinical testing of feasibility and performance[END_REF] .

Conceptually, focal therapy has the potential to minimize treatment-related toxicity without compromising cancer-specific outcome [START_REF] Lazzeri | Focal therapy meets prostate cancer[END_REF] . A current limitation to the proliferation of focal therapy as a viable treatment option is the inability to stage or grade the cancer accurately because of suboptimal imaging capabilities. In the near future, there is no doubt that imaging will detect better cancer foci, especially those with high probability to develop [START_REF] Villers | Dynamic contrast enhanced, pelvic phased array magnetic resonance imaging of localized prostate cancer for predicting tumor volume: correlation with radical prostatectomy findings[END_REF] . As shown in recent studies, multiparametric MRI makes it possible to predict more accurately the location of tumor foci particularly those with high grade [START_REF] Sciarra | Advances in magnetic resonance imaging: how they are changing the management of prostate cancer[END_REF] .

We believe that MR imaging could improve the efficiency of focal therapy and obviously prostate biopsies.

MRI guidance is a challenging task due to the need of MR compatible instruments and the difficulty to process images in real time. In contrast, MRI-TRUS fusion has the advantage of combining pre-operative MR images and intra-operative ultrasound images in a real-time display [START_REF] Ukimura | 3-dimensional elastic registration system of prostate biopsy location by real-time 3-dimensional transrectal ultrasound guidance with magnetic resonance/transrectal ultrasound image fusion[END_REF][START_REF] Hadaschik | A novel stereotactic prostate biopsy system integrating preinterventional magnetic resonance imaging and live ultrasound fusion[END_REF] .

The robot described herein is an attempt at a multimodality robot that could combine the accuracy of robotic needle guidance, precise tracking of prostate motions and deformations as well as MRI-TRUS fusion for facilitated target detection and intervention planning. Indeed, PROSPER is coupled to a MRI-TRUS fusion algorithm. A 3D US reference volume is acquired then registered to a pre-operative MRI acquisition to facilitate and improve prostate delineation, initial dose planning as well as the identification of suspicious areas. Note that MRI fusion was not used in the experiments reported in this paper, because it could interfere with the evaluation of the prostate tracking, which was our main objective. The fusion technique has, however been detailed in previous publications by our group [START_REF] Reynier | MRI/TRUS data fusion for prostate brachytherapy. Preliminary results[END_REF][START_REF] Daanen | MRI/TRUS data fusion for brachytherapy[END_REF] .

As previously mentioned, this system has been described initially for a brachytherapy application. Indeed brachytherapy for low-risk localized prostate cancer is a standard of care, while results for focal therapy require more maturity. With this in mind, a dosimetry planning software taking into account the oblique directions of the needles is in process.

The system could, also, be used for prostate transperineal biopsies, as the described system attempts, in general, to help the surgeon to perform accurate punctures in the operating room with the patient in a lithotomic position. Given that TRUS-MRI fusion is available on our system, it could be useful to perform multiparametric prostate biopsies. However, a transperineal approach would be required to use the robot.

The originality of the PROSPER robot is its ability to track prostate motions and deformations to correct the trajectories and depth of the needles [START_REF] Baumann | Prostate biopsy tracking with deformation estimation[END_REF] . During a manual procedure, especially in brachytherapy, these motions are suspected and the dosimetry plan is manually registered with the prostate shape. This study shows that these motions are not so easily noticeable and that positioning errors may occur also because of prostate rotation, which cannot be evaluated with classic US modalities.

During a prostate brachytherapy procedure, such a robot might be clinically useful to correct for tissue deflection and deformation, adjust seed positions based on real-time dosimetry feedback and avoid the pubic arch so as to extend indications for large prostates. Furthermore, a reduction of the number of needles required in the target volume could be enabled by inverse optimization to create more conformal treatment plans.

These goals can be achieved by oblique needle insertions that require a robotic system. Indeed, manual insertion requires placing needles into areas according to the standard template holes.

Additionally the tracking system has the ability to make fine positional adjustments to offset distribution depending on gland motion.

This robot, by providing a fusion between MRI images and an accurate 3D-TRUS guidance, could represent a future solution for robotic targeted focal therapy [START_REF] Eggener | Focal therapy for prostate cancer: possibilities and limitations[END_REF][START_REF] Reynier | MRI/TRUS data fusion for prostate brachytherapy. Preliminary results[END_REF][START_REF] Long | Prostate biopsies guided by three-dimensional real-time (4-D) transrectal ultrasonography on a phantom: comparative study versus two-dimensional transrectal ultrasound-guided biopsies[END_REF] . In the mean time, the system can be useful for performing more accurate robotic brachytherapy and prostate biopsies.

The study has limitations. It can be considered as a first experience with this new robot. In the statistical analysis of our phantom experiments, it is important to keep in mind that the number of phantoms used was chosen empirically.

In our model, we did not consider the potential bleeding due to the needle insertion in vivo that can usually cause imaging artifacts, however, each needle inserted into the phantom left a prominent trace in the image, progressively degrading US image quality. This degradation did not affect the US registration algorithm, so we are confident that bleeding will not significantly affect system accuracy.

An important future step for the furthered success of our prostate needle insertion system would be to determine the degree to which the prostate rotation occurs in vivo and to provide ways of mitigating this error, such as predicting the motion with biomechanical models or reducing prostate rotation with stabilizing needles [START_REF] Podder | Methods for prostate stabilization during transperineal LDR brachytherapy[END_REF] .

Before using this system on human beings, preclinical studies on more realistic models are needed. These experiments encourage us to continue on with an in depth pre-clinical cadaver study.

CONCLUSION

This robot for prostate focal therapy, brachytherapy and targeted prostate biopsies is the first system to use intraoperative prostate motion tracking to guide needles into the prostate. The preliminary experiments described show its ability to reach targets in spite of the motion of the prostate.

Figure 1 :Figure 2 :

 12 Figure 1: Schematic representation: a-the PROSPER robot design, b-System setup

Figure 3 :

 3 Figure 3: Glass bead implantation in a mobile and deformable phantom a-US view of the beads embedded inside the phantom with definition of a target b-CT scan showing the target hit by a glass bead

  to simulate pubic arch avoidance. The phantoms were then imaged in a Philips™ Brilliance 64 CT scanner at a scanning resolution of 0.15 x 0.15 mm per pixel and 0.33 mm slice spacing. The target and

inserted beads were segmented and the distance between them was measured. Data were presented as median and interquartile range. Phantom coordinates were categorized into locations (anterior-posterior, right-left-central zone). Comparison of the accuracy according to the locations was performed using a Mann-Whitney U test or a Cruskal-Walis test as appropriate with PASW v18 (SPSS Inc, Chicago,Il).

Table 1 : Distance errors between targets and beads stratified according to spatial location # beads Distance error p Depth correction p Apex

 1 Values are expressed as median (IQR). Unit is mm.

50 2.28 (1.8-2.7) <0.001* 4.0 (3.0-5.6) <0.001* Base 40 3.83 (2.9-4.7) 6.5 (5.8-8.2) Center 28 2.6 (2.1-3.4) 0.002* 5.9 (3.1-6.9) 0.6 Left 32 3.6 (2.7-4.8) 7.1 (3.4-11.2) Right 30 2.2 (1.8-2.9) 5.7 (3.8-9.5) Anterior 52 2.7 (2.1-3.6) 0.59 5.4 (3.1-6.7) 0.18 Posterior 38 2.7 (2.1-3.9) 5.9(4.6-6.8) Horizontal 67 2.7 (2.1-3.8) 0.28 5.4 (3.1-6.7) 0.18 Angled 23 2.82 (1.8-3.3) 5.9(4.6-6.8)

Table 2 : Median distance of target movements along each axis showing mobility in all three directions and suggesting, not just prostate translation in the needle direction (z), but also prostate rotation.

 2 Results are expressed as median and IQR (mm)

		x	y	z
		1.26 (0.57-	1.09 (0.51-	1.53 (0.89-
	overall	2.34)	1.78)	2.25)
		1.10 (0.54-	0.73 (0.37-
	Apex	1.81)	1.47)	1.36 (0.7-1.89)
	Base	2.1 (0.81-3.78) 1.49 (0.7-2.14) 2.0 (1.05-2.57)
		0.97 (0.37-	0.89 (0.49-
	Right	1.31)	1.53)	(0.7-2.0)
			1.35 (0.55-
	Center	2.03 (1.20-3.7)	2.17)	1.59 (1.0-2.6)
		2.13 (1.23-
	Left	3.77)	1.45 (0.57-2.3) 1.69 (1.03-2.4)
	Anterior 1.1 (0.55-1.81) 0.73 (0.37-1.7) 1.36 (0.7-1.89)
	Posterio			1.21 (0.69-
	r	1.62 (0.61-2.7) 1.09 (0.62-1.7)	2.11)

This work was funded by the French Agence Nationale de la Recherche (ANR) under the "Technologies pour la Santé" program (Prosper project, coordinated by Koelis SAS, La Tronche, France). J.Troccaz is supported by INSERM (Contrat Hospitalier de Recherche Translationnelle 2010).


Pierre Gosselet

Christian Rey

Daniel J Rixen
email: d.j.rixen@wbmt.tudelft.nl

On the initial estimate of interface forces in FETI methods

Keywords: domain decomposition, iterative solver, FETI, Schur complement, force splitting

The Balanced Domain Decomposition (BDD) method and the Finite Element Tearing and Interconnecting (FETI) method are two commonly used non-overlapping domain decomposition methods. Due to strong theoretical and numerical similarities, these two methods are generally considered as being equivalently efficient. However, for some particular cases, such as for structures with strong heterogeneities, FETI requires a large number of iterations to compute the solution compared to BDD. In this paper, the origin of the bad efficiency of FETI in these particular cases is traced back to poor initial estimates of the interface stresses. To improve the estimation of interface forces a novel strategy for splitting interface forces between neighboring substructures is proposed. The additional computational cost incurred is not significant. This yields a new initialization for the FETI method and restores numerical efficiency which makes FETI comparable to BDD even for problems where FETI was performing poorly. Various simple test problems are presented to discuss the efficiency of the proposed strategy and to illustrate the so-obtained numerical equivalence between the BDD and FETI solvers.

Introduction

Domain Decomposition methods provide a natural framework to solve engineering problems decomposed into subparts. Such problems can arise for instance because each subdomain is discretized independently, or because the subdomains represent different physical domains. Decomposed domains can also be created from an initial single domain problem in order to make efficient use of parallel computing hardware, i.e. in order to distribute the computing work on several processors. Among the domain decomposition methods applied in engineering mechanics to solve elliptic linear problems, two similar methods have emerged in the last decade as efficient parallel computing methods: the primal and the dual Schur complement methods. More specifically, two procedures have been shown to ensure scalability and robustness: the Balanced Domain Decomposition (BDD) method [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Tallec | Domain-decomposition methods in computational mechanics[END_REF] and the Finite Element Tearing and Interconnecting (FETI) method [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. The BDD method is a primal procedure where preconditioned conjugate gradient iterations are applied to find the interface displacements that satisfy the interface equilibrium. In the FETI approach, it is the interface forces that are searched for iteratively so as to satisfy the interface displacement compatibility. Therefore FETI is sometimes referred to as a dual method.

The primal and dual Schur complement methods are based on very similar concepts (see e.g. [START_REF] Rixen | Encyclopedia of Vibration[END_REF] for a mechanical description). Mathematically, it has been shown that the preconditioned interface operators for both the BDD and FETI methods have a condition number bounded by [START_REF] Mandel | Convergence of a substructuring method with Lagrange multipliers[END_REF][START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: Connections and new results[END_REF]

κ = O 1 + log H h 2 (1)
where H and h represent the subdomain and the mesh size respectively. Hence it is often accepted in the Domain Decomposition community that using one method instead of the other is a matter of taste and implementation preferences. However, when looking closely at the details and variants of both methods, it becomes rapidly clear that showing the exact equivalence between the primal BDD and the dual FETI is not easy, if at all possible [START_REF] Fragakis | A unified framework for formulating domain decomposition methods in structural mechanics[END_REF].

Let us then consider the simple example depicted in Figure 1 of a highly heterogeneous elastic cube ( E1 E2 = 10 5 ) subdivided into 3 × 3 × 3 subdomains. A uniform pressure is applied on the face opposite to the clamped side. The structure is discretized using Q 2 hexahedral finite elements with 27 nodes. The model contains 21000 degrees of freedom, of which 6000 belong to subdomain interfaces. In Figure 2 the convergence curves of the global equilibrium residual corresponding to iterations of the BDD and the FETI methods are plotted. Both methods are equipped with what literature refers to as the best preconditioners and coarse grids (see section 2). Although the asymptotic convergence of both procedures is similar, it is observed that the convergence of the FETI method is less monotonic and that its initial residual is significantly higher. These results seem to indicate that some specific details such as the choice of the initial estimates result in possibly significant differences between the primal and dual algorithms.

E1

In this paper we revisit the way the interface forces are estimated initially in the FETI method in an attempt to obtain a convergence comparable to that of the BDD.

In the next section, we shortly recall the concepts underlying the FETI solver.

In section 3 we explain that the forces applied on the interface can be split in different ways. Although the final result is independent of that splitting, it affects the actual FETI iteration history. We then present an efficient way to define such a splitting and show how it is related to the construction of the initial iterate of FETI in section 4. Numerical examples are reported in section 5 to illustrate the effectiveness of the new initialization strategy. Finally, we present some conclusions.

FETI basics

The decomposed problem

Let us consider a domain Ω subdivided into N s non-overlapping subdomains Ω (s) and assume that we are solving a linear (or linearized) static equilibrium problem on the domain. The discretized subdomain equilibrium is expressed by

K (s) u (s) = f (s) + g (s) s = 1, . . . N s (2)
where K (s) , u (s) and f (s) are the subdomain stiffness matrices, displacements and applied forces respectively. g (s) are the connecting forces on the interface between subdomains (thus zero on the internal degrees of freedom). For the sake of simplicity, we assume in the following that the meshes are matching (conforming) on the interface.

The interface forces satisfy an interface equilibrium equation expressing that when assembled on the interface, the resultant is null (action-reaction):

Ns s=1 L (s) T g (s) = 0 (3)
where L (s) is a Boolean assembly matrix. The interface connecting forces are such that the interface degrees of freedom are compatible, namely

Ns s=1 B (s) u (s) = 0 (4)
This relation expresses that for any pair (u (s) , u (r) ) of degrees of freedom matching and the interface, u (s)u (r) = 0. B (s) are thus signed Boolean matrices expressing the compatibility constraints on the interface.

The equilibrium problem of domain Ω is fully described by the local equilibrium (2) and by the interface constraints [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. In block diagonal notations, it can be summarized as

             Ku = f + g L T g = 0 Bu = 0 ( 5
)
where K is the block diagonal matrix of the local operators K (s) and where

u =       
u (1) . . .

u (Ns)        f =       
f (1) . . .

f (Ns)        g =       
g (1) . . .

g (Ns)        L T = L (1) T • • • L (Ns) T B = B (1) • • • B (Ns) (6)
Note that in this description, one set of interface displacements and one set of interface forces are defined per subdomain.

Solvers for decomposed problems

Solving (5) can be done in several ways:

• Considering (5) as a constrained equilibrium problem in terms of u and g leads to the three-field formulation of decomposed domains (see e.g. [START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: Algorithm description[END_REF][START_REF] Rixen | Theoretical comparison of the feti and algebraically partitioned feti methods, and performance comparisons with a direct sparse solver[END_REF]). • One can choose to work with a displacement set u that satisfies a priori the interface compatibility [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. For that purpose we define a global set u g of degrees of freedom unique on the interface such that

u (s) = L (s) u g or u = Lu g (7)
where L (s) is the same assembly Boolean matrix as in (3) that extracts subdomain degrees of freedom from the global set. Stating that u (s) are obtained from a unique set is obviously equivalent to stating the interface compatibility ( 4) and ( 7) thus implies

Bu = BLu g = 0 (8)
for any global displacement u g . On the other hand, all compatible displacements can be written as in [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: Connections and new results[END_REF]. Hence

L = null(B) (9)
In order to illustrate these concepts, we reproduce the example of [START_REF] Rixen | Extended preconditioners for FETI method applied to constrained problems[END_REF] in Figure 3. Introducing [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: Connections and new results[END_REF] in [START_REF] Rixen | Encyclopedia of Vibration[END_REF] yields

     KLu g = f + g L T g = 0 (10)
This set of equations is at the basis of the primal iterative solution techniques such as the Primal Schur Complement or the BDD methods [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Tallec | Domain-decomposition methods for large linearly elliptic three dimensional problems[END_REF]: iteration schemes are applied to find the displacements u g until the interface equilibrium L T g = L T (KLu gf ) = 0 is satisfied. • One can choose in (5) a set of interface forces satisfying a priori the interface equilibrium L T g = 0 while keeping redundant interface degrees of freedom in u. According to [START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: Algorithm description[END_REF], such interface forces have the generic expression

g (s) = -B (s) T λ or g = -B T λ (11)
λ are interface forces that act in opposite directions between any pair of matching degrees of freedom on the interface and are therefore in equilibrium

Ω (1)
Ω ( 2)

Ω (3) u (1) 1 u (2) 1 u (1) 2 u (2) 2 u (3) 1 λ 1 λ 2 λ 3 λ 4 B (1)         -1 0 0 -1 0 -1 0 0 B (2) 1 0 0 1 0 0 0 -1 B (3) 0 0 1 1         B            u (1) 1 u (1) 2 u (2) 1 u (2) 2 u (3) 1            = 0            u (1) 1 u (1) 2 u (2) 1 u (2) 2 u (3) 1 
          =            1 0 0 1 1 0 0 1 0 1            L   u g 1 u g 2  
Fig. 3. Lagrange multipliers and interface compatibility (see Figure 3). Problem (5) becomes

     Ku + B T λ = f Bu = 0 (12)
Clearly, λ are the Lagrange multipliers associated to the interface compatibility constraints. This form of the decomposed problem is the basis for the dual procedures such as FETI: iterative algorithms are applied to compute the interface forces λ such that the displacements resulting from the subdomain equilibrium are compatible on the interface. • If one chooses interface displacements that are unique on part of the interface while, on the remainder of the interface, equilibrated connecting forces are defined, one obtains hybrid primal/dual approaches such as the FETI-DP procedure [START_REF] Farhat | FETI-DP: a dualprimal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF]. • If on the entire interface we use displacements and forces that satisfy a linear combination of the interface equilibrium and compatibility (e.g. Robin type of boundary conditions), one obtains formulations typically used in wave propagation analysis such as described for instance by Helmholtz equations [START_REF] De La Bourdonnaye | Advances in Computational Mechanics with High Performance Computing[END_REF].

• Finally if both the interface equilibrium and compatibility are enforced a priori, one obtains the fully assembled form

L T KLu g = L T f (13)

FETI: the dual iterative solver

In the FETI method [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF], the decomposed problem ( 12) is expressed in terms of interface forces λ: using the subdomain equilibrium equations to eliminate u (s) ,

u (s) = K (s) + f (s) -B (s) T λ -R (s) α (s) (14)
where K (s) + is the inverse of K (s) or a generalized inverse if subdomain Ω (s) is floating when disconnected from its neighbors. In the latter case, R (s) are the associated rigid body modes, their amplitudes α (s) being determined such that the interface forces are in equilibrium with the applied forces f (s) , i.e. such that the subdomain equilibrium is well-posed:

R (s) T f (s) -B (s) T λ = 0 (15)
Substituting ( 14) into the interface compatibility condition, and taking account of ( 15), one obtains the dual interface problem

   F I G I G T I 0       λ α    =    d e    (16)
where

F I = Ns s=1 B (s) K (s) + B (s) T = BK + B T d = Ns s=1 B (s) K (s) + f (s) = BK + f G I = B (1) R (1) • • • B (Ns) R (Ns) = BR α =       
α (1) . . .

α (Ns)        and e =       
R (1) T f (1) . . .

R (Ns) T f (Ns)        = R T f
Let use define the projection operator P such that G T I P = 0:

P (Q) = I -QG I G T I QG I -1 G T I ( 17
)
The choice of the operator Q is discussed later. Introducing the splitting

λ = P (Q) λ + λ 0 ( 18
)
λ 0 = QG I G T I QG I -1 e (19)
the dual interface problem ( 16) is equivalent to

P (Q) T F I P (Q) λ = P (Q) T (d -F I λ 0 ) (20)
The FETI method consists in preconditioned conjugate gradient iterations on the dual interface problem (20) [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. Applying F I to an iterate is a naturally parallel operation. The projection steps however require solving a coarse global problem. The preconditioning operators and the possible choices for Q are summarized next.

FETI preconditioners and coarse grid space

The forces B T λ exist only on the interface degrees of freedom. Using a subscript b and i for interface boundary and internal degrees of freedom respectively, one can re-write the decomposed problem (12) by condensing out u i :

     Su b + B T b λ = f * b B b u b = 0 (21)
where B b is a submatrix of B associated to boundary displacements only and

S =        . . . K (s) bb -K (s) bi K (s) -1 ii K (s) ib . . .        (22) f * b =        . . . f (s) b -K (s) bi K (s) -1 ii f (s) i . . .        (23)
S is the block diagonal matrix of the local operators statically condensed on the interface (also known as Schur complements). The dual interface problem [START_REF] Bhardwaj | Application of the FETI method to ASCI problems: Scalability results on a thousand-processor and discussion of highly heterogeneous problems[END_REF] or [START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF] can thus be written in equivalent forms by replacing K, f and B by S, f * b and B b respectively.

The preconditioner in FETI for which the optimal conditioning number (1) holds is the Dirichlet preconditioner:

F -1 I D = Bb S Bb T (24)
The operator Bb is similar to B b but includes a scaling with respect to interface multiplicity or relative interface stiffness [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF]. The scaling is implemented as a simple pre-and post-processing in the preconditioning step. Mathematically, its general expression is [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: Connections and new results[END_REF][START_REF] Rixen | Extended preconditioners for FETI method applied to constrained problems[END_REF]

Bb = B b AB T b + B b A (25)
where the + superscript denotes an inverse or a pseudo-inverse in case redundant compatibility constraints are present. It can be shown that the multiplicity scaling procedure corresponds to A = I, whereas the stiffness scaling (also known as super-lumped scaling) corresponds to A = diag(K bb ) -1 .

The operator Q in the coarse grid projector [START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF] defines the interface forces QG I associated to rigid mode displacements of subdomains. It can be chosen as Q = I, but for many engineering problems it should be taken as the preconditioner, namely

Q = F -1 I D
(26) or at least as a degenerated form of the preconditioner such as

Q = B b diag(K bb ) -1 B T b + (27)
Further discussion on preconditioning, scaling and on the choice of Q can be found in [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Rixen | Extended preconditioners for FETI method applied to constrained problems[END_REF][START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF][START_REF] Bhardwaj | Application of the FETI method to ASCI problems: Scalability results on a thousand-processor and discussion of highly heterogeneous problems[END_REF].

Splitting forces on an interface

The definition of the decomposed hybrid problem [START_REF] Tallec | Domain-decomposition methods for large linearly elliptic three dimensional problems[END_REF] is not unique: any interface force satisfying the interface equilibrium can be added to the system without changing the final solution u. Indeed, replacing [START_REF] Tallec | Domain-decomposition methods for large linearly elliptic three dimensional problems[END_REF] by

     Ku + B T λ = f = f + B T µ Bu = 0 (28)
for any µ will yield the same solution u, the Lagrange multiplier being then such that λ = λ+µ. This can also be observed from the fact that ( 12) and ( 28) have the same assembled forces since L T f = L T f . Mechanically speaking, it means that adding forces on one side of the interface boundary and subtracting it on the other side yields the same solution in terms of displacement although it modifies the internal forces on the interface.

Force splitting commonly used

From the discussion above, it is clear that one can choose the splitting of interface forces on the interface boundaries1 . Commonly, interface forces are split between the connecting subdomains proportionally to subdomain relative stiffness, in a way consistent with the scaling in the preconditioner. Mathematically speaking, this is equivalent to choosing

f = diag(K)L L T diag(K)L -1 L T f (29)
where L T f are the applied forces given in an initially assembled problem or the assembled forces corresponding to a different force splitting. (L T diag(K)L) is the assembled diagonal stiffness on the interface.

Let us note that for any symmetric positive definite matrix A, the following relation holds:

2 AL L T AL -1 L T + B T BA -1 B T + BA -1 = I (30)
From (30), one deduces that the force splitting (29) is equivalent to

f = f -B T B diag(K) -1 B T + B diag(K) -1 f (31)
This last relation shows that the common splitting technique (29) corresponds to adding a particular set of equilibrated interface forces as described by (28).

Splitting statically condensed forces

In section 2.4, we indicated that the hybrid decomposed problem can be set in an equivalent form condensed on the interface (21). Following the same splitting procedure for the force as in (29), one would then construct the problem

     Su b + B T b λ = f * b B b u b = 0 (32)
where

f * b = diag(K bb )L b L T b diag(K bb )L b -1 (L T b f * b ) (33)
and where (L T b f * b ) are the applied forces statically condensed and assembled on the interface. Following a similar discussion as in the previous section, it is straightforward to show that such a splitting is equivalent to

f * b = f * b -B T b B b diag(K bb ) -1 B T b + B b diag(K bb ) -1 f * b (34)
indicating again that the splitting corresponds to adding a particular set of equilibrated interface forces.

Let us observe that, in the non-condensed format, (32) is identical to

Ku + B T λ =    f i f * b + K bi K -1 ii f i    (35) = f -B T B diag(K) -1 B T + B diag(K) -1    0 f * b   
When comparing it with (31), it is clear that splitting the condensed forces is not equivalent to splitting f as commonly done. Although the displacement field obtained is obviously the same, the Lagrange multipliers searched for during the FETI iterations will be different. The question thus arises: what splitting should be used and how does affect the FETI iterations on the dual interface problem? Next section provides a theoretical analysis of the interest of the splitting while the section after provides related numerical assessments.

An improved initial estimate of the Lagrange multipliers

Let us consider again the hybrid decomposed problem [START_REF] Tallec | Domain-decomposition methods for large linearly elliptic three dimensional problems[END_REF]. In order to construct an estimate for the interface Lagrange multipliers, let us assume that the internal degrees of freedom u i satisfy the local equilibrium while u b on the interface boundary have a zero estimate. This is exactly the assumption underlying the initialization of the primal Schur complement iteration schemes. The subdomain equilibrium then writes

B T b λ ≃ f * b ( 36
)
which cannot be exactly satisfied (unless u b = 0 corresponds to the solution).

Hence, let us decompose the statically condensed force f * b into a component satisfying the interface equilibrium (i.e. belonging to the image of B T b ) and a remainder:

B T b λ ≃ f * b = B T b γ + remainder (37)
λ 00 = B b diag(K bb ) -1 B T b + B b diag(K bb ) -1 f * b (40)
can be easily computed by a scaling of the interface condensed forces. Observing that diag(K bb ) -1 is an approximation of S + , λ 0 minimizes a norm having the meaning of an energy.

We thus propose to start the FETI iterations with an initial estimate obtained by rendering λ 00 admissible, i.e. such that G T I λ 0 = e. Using the notations introduced in (17),

λ 0 = P (Q)λ 00 + QG I G T I QG I -1 e (41)
To be consistent with the choice D = diag(K bb ) -1 , one should take Q as in ( 26) or ( 27).

Remarks

•

If D = S + in (38), λ 00 = F I + d. Furthermore, one would have Q = F I + and λ 0 = P (Q)F I + d + F I + G I G T I F I + G I -1 e = F I + (d -G I α) (42) with α = G T I F I + G I -1 (G T I F I + d -e)
Hence the new choice of initial estimate would yield the exact solution. • If we choose to split the interface forces such as described in (33) (or its equivalent form 34), (40) yields λ 00 = 0 so that the initial iterate (41) would correspond to the standard FETI starting procedure. Hence it is equivalent to apply standard FETI iterations when splitting condensed forces on the interface as in (33) or to use any decomposed force vector f together with the starting procedure (41). • The construction of λ 00 in (40) has exactly the same mechanical interpretation as the scaling of interface forces in the preconditioner (see [START_REF] Rixen | Extended preconditioners for FETI method applied to constrained problems[END_REF]). The cost incurred by ( 33) or (41), corresponds to a preconditioning step, thus less than one half of a FETI iteration. Obviously, if the lumped preconditioner is applied, the standard splitting (29) should be considered.

Numerical assessment

In order to assess the performance of the new estimate for the Lagrange multipliers, we consider the solution to various problems by FETI and BDD. In all cases, convergence is monitored through the evaluation of the global residual:

K g u g -f g f g ε (43)
where K g , u g and f g are global assembled stiffness matrix, displacement field and forces (see [START_REF] Farhat | FETI-DP: a dualprimal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF]). All results are obtained with the Dirichlet preconditioner equipped with the super-lumped scaling, i.e. A = diag(K bb ) -1 . The stopping criterion ε is set to 10 -6 . The FETI method is tested for its different projectors: we note P ( F -1 I D ) the Dirichlet projector, P (W ) the superlumped projector corresponding to (27) and P (I) the identity projector. 5.0.0.1 Cube with checkerboard heterogeneity: Let us first assess the strategies described above on the problem described in the introduction (see fig. 1). Figure 4 presents the convergence of the primal residual through the conjugate gradient iterations for the BDD and for the FETI approaches with standard force splitting. The primal approach clearly exhibit better results. Also, due to the strong heterogeneity in the structure and along the interfaces, the projector P (I) yields very poor convergence and is thus not suitable. It is observed that, although the Dirichlet projector yields a better convergence rate, the superlumped projector P (W ) has a significantly lower initial residual and hence reaches convergence faster. The initialization associated to P (W ) leads to a starting residual 10 3 times smaller, however its behaviour during the solution process is less regular. Figure 5 presents the convergence history for the BDD method and for the FETI method equipped with the Dirichlet projector and with the new initialization. For comparison, the convergence obtained with the standard splitting is shown once more. Note that the new initialization not only provides a better initial estimate (lower initial residual) but also leads to similar convergence rates as for the standard splitting. Therefore the convergence history of FETI becomes very similar to the convergence of the primal BDD approach. The total number of iterations is about the same as for the primal BDD, although the convergence of FETI is less monotonic. 1 Performance results for the cube of fig. 1 Table 1 summarizes the performance results of the various available strategies. To investigate the efficiency of the new initialization (or force splitting) procedure, we compare the number of iterations to achieve convergence and the norm of the initial primal residual. The new initialization yields significant improvements for the Dirichlet projector (35% less iterations), but only slightly affects the convergence when other projectors are used. 5.0.0.2 Cube and slanted cube with heterogeneous layers: We now assess the new initialization strategy on two other configurations of the cube in order to evaluate the influence of the geometry and of the repartition of heterogeneities. The structures depicted in Figure 6 and 7 are similar to fig. 1 but with different material distribution. Also, for the problem described in Figure 7, the cube has been slanted by 60 degrees. In table 2 we report the number of iterations when using the BDD solver and when applying FETI with the standard and the new initialization. For the straight cube of fig. 6, the use of the novel initialization leads to only small improvements. Observing that the FETI solver with Dirichlet projector and standard splitting converges nearly as fast as the BDD, it is clear that the new splitting strategy can not have a significant effect.

For the slanted cube, due to the geometric distortion of the substructures and of the mesh, the BDD convergence is better than when the FETI with the standard splitting is applied. When the new initialization strategy is used, the convergence history of FETI is again very similar to the convergence of the primal BDD. 5.0.0.3 Nonlinear flexion of a composite beam As a last example, we analyze a non-linear problem solved through a sequence of linearized systems. The structure (fig. 8) is a slender beam of aspect ratio 9 with square crosssection. It is made of longitudinal strips of metal and rubber. The loading corresponds to an imposed pressure on one side of the beam. Due to the presence of elastomer parts, the structure undergoes large deformations. We chose 6 and7 a Kirchoff Saint-Venant model for the metal, the characteristic coefficients of which are its Young modulus E = 20000 MPa and Poisson's coefficient ν = 0.3. A NeoHookian model is assumed for the rubber characterized by a shear modulus G = 2.0 MPa and a compressibility modulus K = 2000 MPa. The structure is decomposed into 27 monomaterial paralellepipedic subdomains and each substructure is meshed in 2 × 2 × 18 cubic elements. In order to handle the quasi-incompressibility of the elastomer, a mixed finite element formulation is considered where the pressure field and displacements are discretized independently. We choose a Q 2 -P 1 hexaedral element with 27 displacement nodes and 4 internal pressure nodes. Details on this formulation and on the practicalities for applying it properly in simulation can be found in [START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF].

The complete model has 55300 degrees of freedom of which 16400 belong to the interface. Due to the behaviour of the elastomer, the problem is highly nonlinear. Newton-Raphson iterations are performed where the tangent matrix is updated at every step. The pressure loading is 5 bars. The stopping criterion is set to 10 -3 for the relative primal residual of FETI when solving the linearized systems. For the Newton-Raphson iterations, the tolerance for the relative residual of the nonlinear equations is set to 10 -4 , so that 5 Newton-Raphson iterations and thus 5 linear solves must be performed. and the Dirichlet preconditioner with stiffness scaling is used. The significant increase of the number of iterations between the first and the following linearized systems can be explained by the loss of positivity of the tangent matrix which is mostly due to incompressibility. Indeed, the conjugate gradient algorithm (with full re-orthoganilization) applied on the interface problem in the FETI and in the BDD methods remains applicable for non positive matrices but its convergence is significantly slowed down [START_REF] Paige | Approximate solutions and eigenvalue bounds from krylov subspaces[END_REF]. 3 As observed from figure 9, the new initialization enables FETI to achieve performance results which are very similar to the BDD method. Comparing the total number of iterations, the classical FETI approach requires 10% more iterations than BDD while the FETI method with new initialization requires slightly less iterations than the BDD. From figure 9 we also observe that the new initialization technique improves the performance of FETI both for the Dirichlet and the diagonal projector. For this particular structure with regular geometry, the cost effective diagonal projector yields a convergence rate very similar to the convergence rate obtained with the more computationally intensive Dirichlet projector, except for the very first linearized system solve.

Our numerical experiments indicate that for a large class of nonlinear problems such as the one depicted here, the new initialization leads to a small but non-negligible gain in terms of number of iterations and CPU time. Another important beneficial effect of the new starting procedure for FETI comes from the fact that, since the initial residual of the iterations on the interface problem are several orders of magnitude lower with the new initialization, stagnation of the residual of the FETI iterations which often happens when dealing with higher nonlinearity (higher loading) is significantly delayed, so that in practice restarting of the iteration can be avoided. For some particular structural problems such as those exhibiting strong heterogeneities and geometric distortion, the Finite Element Tearing and Interconnecting (FETI) solver can yield poor convergence compared to the conceptually similar Balanced Domain Decomposition (BDD) solver. In those cases, the bad performance of FETI can be traced back to high initial residual in the iterations on the interface problem. The initial residual is strongly related to the way the applied forces are split on the subdomain interface boundaries.

In this paper, we propose a novel strategy to split the applied forces between the subdomains. We propose to split the statically condensed interface force according to interface diagonal stiffness. This leads to building a more efficient initial estimate for the interface connecting forces. The new initialization for the FETI iterations mainly involves computing statically condensed forces on the interface and can thus be performed at a computational cost equivalent to less than half the cost of a full FETI iteration.

The numerical examples described in this paper indicate that for the problems where the primal BDD method outperforms FETI due to unexpected high initial residual, the new initialization strategy builds a better starting estimate of the interface forces and, in turn, to an initial residual similar to the BDD residual. The FETI method then converges in a manner very similar to the BDD method.

The novel starting strategy never deteriorates the FETI convergence and leads to significant improvements in some pathological cases. Therefore we suggest to use the presented initialization as default in FETI solvers.

With the proposed initialization for FETI, the FETI method and the BDD solver lead to similar convergence and computational costs for complex problem where the Dirichlet projector is required. For problems where the simplified FETI preconditioners and projectors can be used without significantly deteriorating the convergence of the interface iterations, FETI is often found to be more efficient in terms of overall computing cost.

Fig. 2 .

 2 Fig. 1. Decomposed heterogeneous cube

  -singular and symmetric matrix D. The remainder is then orthogonal to DB T b . If we choose as initial estimate λ 00 = γ (39) then the initial equilibrium residual B T b λ 00f * b = remainder is minimum in the D-norm. Setting D to be diag(K bb ) -1

Fig. 4 .

 4 Fig. 4. Cube with checkerboard heterogeneity (fig. 1): convergence of the primal residual for FETI with standard force splitting

Fig. 5 .

 5 Fig. 5. Cube with checkerboard heterogeneity (fig. 1): convergence of the primal residual for FETI with new initialization

ig. 7 .

 7 Fig. 6. Cube with heterogeneous layers

Figure 9 Fig. 8 . 3 × 3 ×

 9833 Figure9indicates the number of iterations required when solving the linearized systems by the classical FETI method, by FETI with improved initialization and by the primal BDD. The Dirichlet and diagonal projectors are applied,

Fig. 9 .

 9 Fig. 9. Performance results for the flexion of the composite beam 6 Conclusions

In the same manner, it was noted in[START_REF] Rixen | Extended preconditioners for FETI method applied to constrained problems[END_REF] that coefficients of additional constraints can also be arbitrarily split on the interface, which led to the construction of efficient preconditioners for problems with multipoint constraints.

This relation expresses the complementarity of the primal and dual scaling and was already suggested in[START_REF] Fragakis | A unified framework for formulating domain decomposition methods in structural mechanics[END_REF]. Its proof can be obtained by noting that L T (30) and BA -1 (30) are trivially satisfied. Owing to the fact that A is symmetric positive definite and because the image of L corresponds to the nullspace of B, the final result is obtained.

Often, non-positivity is due to the mutation of former null-modes (rotations) to negative modes, it can be handled by the introduction of these modes as constraints in Krylov-augmented algorithms[START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF]. In our case, non-positivity is mostly due to the behaviour of the rubber and the strategy described above is non-relevant. An efficient strategy based on the approximation of negative modes can be found in[START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF].

Acknowledgements

Part of this work was performed while the third author was visiting the LM2S with the financial support of the Ecole Normal Supérieure de Cachan, France. The first two authors acknowlegde computational resources support from the Centre Informatique National Enseignement Supérieur (CINES) and the Pôle de Calcul Paris Sud.