Nicolas Halbwachs

Julien Henry

When the decreasing sequence fails ⋆

The classical method for program analysis by abstract interpretation consists in computing a increasing sequence with widening, which converges towards a correct solution, then computing a decreasing sequence of correct solutions without widening. It is generally admitted that, when the decreasing sequence reaches a fixpoint, it cannot be improved further. As a consequence, all efforts for improving the precision of an analysis have been devoted to improving the limit of the increasing sequence. In this paper, we propose a method to improve a fixpoint after its computation. The method consists in projecting the solution onto well-chosen components and to start again increasing and decreasing sequences from the result of the projection.

Introduction

Program analysis by abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] consists in computing an upper approximation of the least fixpoint of an abstract semantic function in a suitable abstract lattice of properties. When the abstract lattice is of infinite depth, the standard approach [START_REF] Cousot | Static determination of dynamic properties of programs[END_REF][START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] consists in computing an increasing sequence whose convergence is forced using a widening operator ; then, from the obtained limit of the increasing sequence, one can improve the solution by computing a decreasing sequence, by iterating the function without widening. The decreasing sequence may either stop at a fixpoint of the semantic function, or be infinite, but since all its terms are correct solutions, one can stop the computation after a fixed number of terms, or limit its length using a narrowing operator .

Of course, the precision of the result depends both on the ability of the widening operator to "guess" a precise limit of the increasing sequence, and on the information gathered during the decreasing sequence. Intuitively, the increasing sequence extrapolates the behaviour of the program from the first steps of its execution, while the decreasing sequence gathers information about the end of the execution of the program, its loops, or more generally, the way the strongly connected components of its control flow graph are left.

While significant efforts have been devoted to improving the precision of the limit of the increasing sequence (see §1.2 for a quick survey), little attention

⊥ ⊥ ⊥ [0, 0][0, ∞] [0, 0][0, 99] [0, 0][1, 100] [0, 0][100, ∞] [1, 1][100, ∞] ⊥ 2 3 [0, 1] [0, 1] [0, 1][0, 0] [0, ∞][0, ∞] [0, ∞][0, 99] [0, ∞][1, 100] [0, ∞][100, ∞] [1, ∞][100, ∞] [100, ∞][100, ∞] 4 [0, ∞][0, ∞] [0, 99][0, 0] [0, 99]
[0, ∞] i<100? j<100? Fig. 1. Example 1 -nested loops has been paid to the decreasing sequence. It is generally admitted that, when the decreasing sequence reaches a fixpoint, it cannot be improved. However, it appears that such a fixpoint can be far from the least fixpoint, so improving it may have a significant influence. More specifically, we shall see in §1.1 that slight modifications in a program may have a surprising influence on the amount of information gathered during the decreasing sequence.

Motivating examples

Let's illustrate the problem with two very simple examples:

Example 1 is a classical example of what can be obtained by interval analysis [START_REF] Cousot | Static determination of dynamic properties of programs[END_REF]. Fig. 1.a shows the control-flow graph (CFG) of a very simple loop incrementing a variable i from 0 to 100. The increasing sequence consists of 2 iterations; at iteration 2, the widening is applied, and the sequence converges. Iteration 3 shows the descending sequence, which reaches a fixpoint in 1 step. The results are the best possible: i ∈ [100, 100] at the end of the loop. Now, the CFG shown in Fig. 1.b is obtained from the preceding one by nesting a second loop on another variable j within the first one. The increasing sequence converges after 4 steps. Again, the descending sequence reaches a fixpoint in 1 step, but now, the result for i is imprecise: i ∈ [100, ∞] at the end. The reason is that the nested loop neither modifies nor tests the variable i; so, as soon as its interval has been widened to [0, ∞], it will remain unchanged in the inside loop. Notice that it is also the case if we select both loop heads as widening nodes (see Note 1 in §2 on the selection of widening nodes).

Example 2: Our second example illustrates a situation which occurs commonly in reactive programs, cyclically sampling their environment. A first program (Fig. 2.a) is just an infinite loop with a counter modulo 60. It is properly analysed after 2 steps of increasing sequence and one descending step.

[1, 1] [1, 1] [0, ∞] [0, 0] [1, 60] [0, 60] [0, 60] [0, 0] [1, 60] [0, 60] F T n:=0 n:=0 n<60? n:=n+1 iterations ⊤ [0, 0] [0, 0] ⊥ 1 2 3↓
convergence fixpoint (a) An infinite loop with a counter modulo 60

T T F T [1, 1] [1, 1] [0, ∞] [1, 60] [0, 60] [0, 0] [0, ∞] [0, ∞] F n:=0 read(sec)? n:=0 n<60? n:=n+1 n:=0 n<60? n:=n+1 iterations ⊤ [0, 0] [0, 0] 1 2 ⊥ [0, 1] fixpoint [0, 0]
(b) A loop counting the occurrences of "seconds"

Fig. 2. Example 2 -intermittent counting

Now, assume that we don't want to count all loop iterations, but only when some external event (e.g., a "second") is detected. This is done by the program of Fig. 2.b. As before, the increasing sequence converges after 2 steps, but now, the limit is a fixpoint, so there is no decreasing sequence, and the upper bound of the counter is missed. This second example can be simply explained as follows: let (L, ⊑, ⊔, ⊓, ⊤, ⊥) be the abstract lattice, and F be an abstract semantic function from L to L. Let G = Id ⊔ F (i.e., λX.X ⊔ F (X)). Then, it is easy to see that G has the same least fixpoint as F , but G is extensive (i.e., ∀X ∈ L, X ⊑ G(X)). As a consequence, any postfixpoint of F is a fixpoint of G. So, while the limit of the increasing sequence with F may be a strict postfixpoint of F -which can be improved by a decreasing sequence -, this limit will be a fixpoint of G, meaning that there is no decreasing sequence with G. This is exactly what happens with our example, where the dummy branch in the loop of Fig. 2.b adds an identity term to the semantic function at the widening node.

Related works

Beside researches proposing new abstract domains, many existing works aim at fighting the imprecision of analysis, considered to be essentially due to the widening operation. Apart from proposals of systematic design of widening and narrowing operators [START_REF] Cortesi | Widening and narrowing operators for abstract interpretation[END_REF], one can distinguish at least three big tracks: (1) designing smart widening operators, generally dedicated to some specific domains;

(2) avoiding or minimising the use of widening, by focusing either on some classes of programs or on some classes of abstract domains; (3) applying widening and narrowing using smart strategies.

Smart widening. Several proposals concern smart widening operators, especially for the polyhedra domain [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Bagnara | Precise widening operators for convex polyhedra[END_REF]. Widening up to or with thresholds [Hal93,HPR97,BCC + 03,LCJG11] consists in choosing -generally from the conditions appearing in the program -some tentative limits to the widening. Widening with landmarks follows the same idea, but the selection of limits is made dynamically. Widening with care set [START_REF] Wang | Using counterexamples for improving the precision of reachability computation with polyhedra[END_REF] makes use of a proof objective. Some of these proposals can properly deal with some of our examples, mainly because they can reach a precise solution at the end of the increasing sequence. However, they are not independent of the considered abstract domain. Our method will work for any abstract domain, and is compatible with any widening operation.

Avoiding widening. Other authors try to avoid the use of widening. Acceleration techniques [BGP97,WB98,CJ98,BFLP03,GH06] are dedicated to some classes of programs or loops, the effect of which can be exactly computed. Other approaches can be applied only with some kinds of domains -namely "weakly relational domains" [START_REF] Miné | Weakly relational numerical abstract domains[END_REF] or "templates" [SSM04,SSM05] -in which policy iteration [SW04,CGG + 05,GS07] allows least fixpoints to be precisely computed. These methods generally solve our problem, but they are restricted either to some class of programs or to some abstract domains.

Widening strategies. An obvious way of improving the precision of the widening is to delay its application [Hal93,BCC + 03], i.e., applying it only after a fixed number of exact steps or intermittently, or applying it after some loop unrolling [START_REF] Goubault | Static analyses of the precision of floating-point operations[END_REF][START_REF] Putot | Static analysis-based validation of floating-point computations[END_REF]. Some strategies adapt the application of the widening according to the discovery of new feasible paths [Hal93,HPR97,BCC + 03] in the program. In particular, [START_REF] Gopan | Lookahead widening[END_REF][START_REF] Gopan | Reps. Guided static analysis[END_REF] proposes a very clever strategy, called lookahead widening, where a succession of increasing-decreasing sequences are computed for more and more feasible paths of the program. Stratified analysis [START_REF] Monniaux | Stratified static analysis based on variable dependencies[END_REF] is a succession of analyses concerning more and more variables, according to their dependencies. None of these strategies provides a general solution to our problem.

Anyway, while some of these methods can work on some of our examples, none of them specifically address the problem of improving the result of the decreasing sequence. We don't pretend our method is better than these works, but that it is different and complementary.

Contribution and summary

In this paper, we propose a method starting from the result of the decreasing sequence, and trying to improve it as follows: the solution will be projected on some of its components, whose propagation is likely to provide a more precise solution, according to some criteria. New increasing and decreasing sequences will be started from the result of the projection, providing a new solution which can be intersected with the previous one. The method is independent from the abstract domain. We'll show that it properly solves our running examples, and that in some cases, it can gather non trivial information about the end of program execution.

The paper is organised as follows: Section 2 introduce the necessary definitions and notations; our method is presented in Section 3 and illustrated on an example in Section 4; Section 5 proposes some ways for improving the performances and Section 6 gives some experimental results.

Definitions and notations

Abstract lattice. As said before, we assume that the analysis makes use of an abstract complete lattice (L, ⊑, ⊔, ⊓, ⊤, ⊥). We assume this lattice to be of infinite depth. The lattice operations are supposed to be available, together with a widening operator ∇, and the interpretation of each program statement s as a function (predicate transformer) f s : L → L.

Control-flow graph. A control-flow graph (CFG) is a graph (N, E), where

the finite set N = {ν 1 , ..., ν k } is made of 3 types of nodes: the start nodes, the junction nodes, and the statement nodes. With each statement node ν i is associated a function f i : L → L.

-E ⊆ N ×N is the set of edges. Start nodes have no incoming edge, statement nodes have one incoming edge, junction nodes have several incoming edges.

Remark: for simplicity, each node has a single output. This means that the classical "test nodes" (used in Figures 1 and2) are split into pairs of statement nodes, whose associated function returns an abstraction of the intersection of its argument with the condition of the test ("then" part) or its negation ("else" part). As an example, Fig. 3 shows the CFG of our example 1.b. Semantic equations. The analysis will associate with each node ν i of the CFG an abstract value X i ∈ L, these abstract values being defined by a system of recursive equations:

∀i = 1..k, X i =    ⊤ if ν i is a start node f i (X j ) if ν i is a statement node and (ν j , ν i ) ∈ E (νj ,νi)∈E X j if ν i is a junction node
We'll often write this system of equations as a vectorial fixpoint equation:

X = F (X) in the lattice L k .
Increasing sequence. Since the lattice L is of infinite depth, the Kleene sequence X 0 = ⊥ k , X ℓ+1 = F (X ℓ ) may be infinite. The classical approach consists in computing the increasing sequence

Y 0 = ⊥ k , Y ℓ+1 = Y ℓ ∇F (Y ℓ )
; from the properties of the widening operator, this sequence is guaranteed to converge after a finite number of steps towards a limit Y ∇ , which is a postfixpoint of F , i.e., F (Y ∇ ) ⊑ Y ∇ . Of course, the increasing sequence is computed in a chaotic way (cf. Figures 1 and2), by propagating changes along the paths of the CFG, and since the widening operation loses information, it is only applied on a selected set W of widening nodes intersecting each loop of the CFG.

Note 1 (on the selection of widening nodes). The set W of widening nodes must be as small as possible, to minimise the number of applications of the widening operator. Since finding a minimal cutting set W is an NP-complete problem, the heuristic classically applied is the method of strongly connected subcomponents (SCSC) proposed by Bourdoncle [Bou93]: the method recursively uses Tarjan's algorithm [START_REF] Tarjan | Depth-first search and linear graph algorithms[END_REF] to find the strongly connected components (SCC) of a directed graph, together with an entry node to each SCC. Entry nodes are the target of back edges, so they are all junction nodes. Bourdoncle's method consists in adding all SCC entry nodes to W , then removing them from the graph and recursively apply Tarjan's algorithm to the rest of each SCC. The result is a hierarchy of SCSC, each of which being cut by a junction node in W . An obvious improvement of this method (which we did not find published anywhere) consists in considering again the hierarchy of SCSC bottom-up, checking whether each SCSC is disconnected by the cut-points of its children. For instance, on the CFG of Fig. 3, a first application of Tarjan's algorithm finds one non-trivial SCC, c 1 = {ν 3 , ν 4 , ν 6 , ν 7 , ν 8 , ν 9 , ν 10 , ν 11 }, with entry node ν 3 . Removing node ν 3 and applying again the algorithm provides the SCSC c 2 = {ν 7 , ν 8 , ν 9 }, with entry node ν 7 , whose removal disconnects the graph. So, Bourdoncle's method provides W = {ν 3 , ν 7 }. Now, since the cut-point ν 7 of the leaf SCSC c 2 also disconnects the father SCSC c 1 , it's enough to choose W = {ν 7 }, as done in Fig. 3.

Decreasing sequence. The limit Y ∇ of the in- creasing sequence is a post-fixpoint of F . If it is a strict post-fixpoint (F (Y ∇ ) < Y ∇ ),
it can be improved by computing a decreasing sequence Z 0 = Y ∇ , Z ℓ+1 = F (Z ℓ ). This sequence can be infinite, or reach a fixpoint of F . In practice, it generally reaches a fixpoint after very few steps. Anyway, since all of its terms are post-fixpoints of F , hence correct approximations of the least fixpoint, one can stop it after a fixed number of steps, or force its convergence using a narrowing operator. In the following, we'll note Z ∆ the last term of the descending sequence, and we'll generally assume that Z ∆ is a fixpoint; however, the results still hold if Z ∆ is a strict post-fixpoint.

(X ℓ ) ⊤ PostFix(F ) (Z ℓ ) Y ∇ ⊥ lfp(F ) Z ∆ (Y ℓ ) Fix(F )
The above figure classically illustrates the sequences in the abstract lattice: (Xℓ) is the (generally infinite) Kleene's sequence, (Y ℓ ) is the (finite) increasing sequence, leading to the post-fixpoint Y ∇ , and (Z ℓ ) is the decreasing sequence of post-fixpoints providing a solution Z ∆ .

3 Improving a (post-)fixpoint solution

An intuition of the solution

Let's look again at example 1.b (see Fig. 3). At the widening node ν 7 , during the decreasing sequence, we have Z 7 = Z 6 ⊔Z 9 . At the end of the decreasing sequence, we find

Z ∆ 7 = (i ∈ [0, ∞], j ∈ [0, 100]), while Z ∆ 6 = (i ∈ [0, 99], j ∈ [0, 0]) and Z ∆ 9 = (i ∈ [0, ∞], j ∈ [1, 100]). Obviously, i ∈ [0, 99], found in Z ∆ 6 , is a correct invariant, which is lost in Z ∆
7 because of the least upper bound with i ∈ [0, ∞] imprecisely found in Z ∆ 9 . Our idea is to start again a propagation of Z ∆ 6 after resetting Z ∆ 9 to ⊥.

Generalised sequences

Restarting an iteration from an arbitrary point requires some changes in the definition of the iteration sequences. We must ensure that a widened sequence starting form an arbitrary point (not necessarily a pre-fixpoint) is increasing; moreover, widening operators are generally designed under the assumption that their first operand is smaller than the second one. We introduce the following notations: Let F be a monotone function from L to L. Let X ∈ L. Then, -we note F ∇ (X) the limit of the sequence

Y 0 = X, Y ℓ+1 = Y ℓ ∇(X ⊔ F (Y ℓ )); -we note F ∇∆ (X) the last term Z ∆ of a descending sequence (Z ℓ ) starting at Z 0 = F ∇ (X)
Remarks:

-The second operand X ⊔ F (Y ℓ ) of the widening is always greater than the first one, and the increasing sequence (Y ℓ ) is indeed increasing. Obviously, F ∇ (X) is the classical approximation of the least fixpoint of the function λX.(X ⊔ F (Y )), i.e., of the least fixpoint of F greater than X. -Notice also that, with X = ⊥, these definitions of sequences and limits match the classical ones recalled in §2. -For any X, F ∇∆ (X) is a correct approximation of the least fixpoint of F , i.e., ∀X ∈ L, F ∇∆ (X) ⊒ lfp(F ). -Neither F ∇ nor F ∇∆ is increasing. As a consequence, there can be some X such that F ∇∆ (X) < F ∇∆ (⊥), i.e., such that the limit obtained from X is more precise than the one computed by the classical iteration.

We address the problem of analysing an SCC of the graph, since the analysis of a complex graph considers each SCC in turn. We consider first the case of an SCC with only one widening node, before addressing the general case.

Case of a single widening node

We need some additional definitions: Path transformers. Let ν i , ν j be two nodes. Let P i,j denote the set of nodes belonging to an elementary path in the CFG going from ν i to ν j . Intuitively the path transformer from ν i to ν j is the function F i,j : L → L, which, from an abstract value X associated with ν i , provides the abstract value F i,j (X) corresponding to the propagation of X along the elementary paths from ν i to ν j . We have:

F i,j (X) =              ⊥ if P i,j = ∅ X if ν i = ν j f j (F i,k (X)) if ν j is a statement node and (ν k , ν j ) ∈ E (ν k , νj ) ∈ E ν k ∈ Pi,j F i,k (X) if ν j is a junction node
Now, let us consider an SCC with only one widening node, ν i . ν i is a junction node (from the way widening nodes are selected). The abstract value at ν i depends on those at the preceding nodes in the SCC, and on abstract values propagated from outside (start nodes and preceding SCCs), which we note Y 0 i since it is the first value = ⊥ at node ν i during the increasing sequence. Let ν j1 , . . . , ν jm be the source nodes of incoming back edges to ν i .The semantic equations considered during the descending sequence can be subsumed as:

ν i Y 0 i ν j1 ν jm Z i = Y 0 i ⊔ Z j1 ⊔ . . . ⊔ Z jm , Z j ℓ = F i,j ℓ (Z i ), ℓ = 1..m
At the end of the sequence, we have

Z ∆ i ⊒ Y 0 i ⊔ Z ∆ j1 ⊔ . . . ⊔ Z ∆ jm (an equality if Z ∆ is a fixpoint).
Projection according to improving components. The idea is to select a set S of nodes, such that the propagation of terms {Z ∆ j ℓ | ν j ℓ ∈ S} is likely to improve the solution. More precisely, we define Z ∆ ⇓ S by

Z ∆ ⇓ S k = Z ∆ k if ν k ∈ S ⊥ otherwise
The choice of S should make F ∇∆ (Z ∆ ⇓ S) more precise or incomparable with Z ∆ , so that Z ∆ ⊓ F ∇∆ (Z ∆ ⇓ S) is an improved result.

For instance, and following the intuition of §3.1, in our example 1.b (Fig. 3), we should choose S = {ν 6 } and restart increasing and decreasing sequences from Z ∆ ⇓ S, i.e., a vector U such that

U k = (if k = 6 then ([0, 99], [0, 0]) else if k = 1 then ⊤ else ⊥)
With this choice, F ∇∆ (U ) is the best possible result (the least fixpoint of F ), shown by the opposite vector.

i j ν 3 [0,100] ν 7 [0,99] [0,100] ν 9 [0,99] [1,100] ν 11 [1,100] [100,100] ν 5 [100,100]
The opposite figure shows the new sequences: the classical solution Z ∆ is projected on U = Z ∆ ⇓ S (generally not a pre-fixpoint), from which a new increasing sequence (U ℓ ) and a new decreasing sequence (V ℓ ) provide a new limit V ∆ . The improved solution is the greatest lower bound Z ∆ ⊓ V ∆ .

Choice of improving components. A component Z ∆ j ℓ is likely to improve the solution if it is strictly smaller than Z ∆ i . It may happen for two reasons:

either some "initial states" in Y 0 i have been left on the paths from ν i to ν jℓ ; this case is not interesting since Y 0 i clearly consists of states which won't be shown to be unreachable. They will belong to any solution.

⊤ (Z ℓ ) Y ∇ ⊥ new lfp(F ) solution Z ∆ U ∇ (U ℓ ) (Y ℓ ) (V ℓ ) V ∆ U = Z ∆ ⇓ S
or, during its propagation along these paths, Z ∆ i has been "truncated" by some condition; this is the interesting case which can add some information. As a consequence, ν j ℓ will be selected in S if

Y 0 i ⊔ Z ∆ j ℓ < Z ∆ i ( Criterion 1)
Moreover, it is useless to propagate again Y 0 i , which already provided the existing result Z ∆ . So, ν j ℓ will be selected in S only if

Z ∆ j ℓ ⊑ Y 0 i (Criterion 2)
Now, let's consider our example 2.b. At convergence, no predecessor of the widening node satisfy our criteria. However, the widening node is preceded by a succession of two junction nodes, the first one being associated with an obviously interesting invariant: n ∈ [0, 60]. A simple change in the CFG taking into account the associativity of junction, would bring this node in the predecessors of the widening node. So, it can be useful to look for "improving nodes" further upstream the widening node. We change our selection process as follows:

A node ν j will be selected in S if it precedes a junction node (possibly

ν i ) (C0) -Y 0 i ⊔ F j,i (Z ∆ j ) < Z ∆ i (C1) -F j,i (Z ∆ j ) ⊑ Y 0 i (C2)
The criterion (C0) above comes from the fact that, during the descending sequence, only junction nodes lose information. The other two criteria generalise our preceding Criteria 1 and 2, by allowing the candidate node ν j not to be an immediate predecessor of ν i . Note that, in our example 2.b, the selected node ν j is such that F j,i = Id. However, our criteria allow also some statement nodes to be on the path from ν j to ν i .

General case

The case of an SCC with several widening nodes is very similar. Only the definition of P i,j needs to be modified: it denotes the set of nodes belonging to an elementary path in the CFG going from i to j without going through a widening node: only the extremities ν i and/or ν j can belong to W . Our criteria for selecting the improving nodes are essentially unchanged: a node ν j will be selected in S if it precedes a junction node (C0) -there exists a widening node ν i such that

• Y 0 i ⊔ F j,i (Z ∆ j ) < Z ∆ i (C1) • F j,i (Z ∆ j ) ⊑ Y 0 i (C2
) and the definition of Z ∆ ⇓ S is unchanged.

A more illustrative example

Our two running examples are properly analysed with the proposed method. They would be also solved with a smart choice of "thresholds" [Hal93,BCC + 03,LCJG11] for limiting the widening during the increasing sequence. Let us consider now another example to show that our method can discover constraints that don't appear as conditions in the program. The example is a rather ad-hoc modification of Example 1.b by variable change.

i := 0; while i < 4 do { j := 0; while j < 4 do { i := i + 1; j := j + 1;} i := ij + 1; Fig. 4 shows the corresponding CFG together with the set of variable states traversed by the execution at ν 7 . Fig. 5 shows the abstract values at widening node ν 7 , during a polyhedra analysis:

} i := 0 i < 4? i ≥ 4? j := 0 j ≥ 4? j < 4? i := i -j + 1 j := j + 1 i := i +
at the end of the classical increasing sequence, we get (0 ≤ j ≤ i); -the decreasing sequence reaches a fixpoint Z ∆ in one step, giving (0 ≤ j ≤ i, j ≤ 4); -at node ν 6 , we have Z ∆ 6 = (0 ≤ i ≤ 3, j = 0), which satisfies all our criteria; we start again from Z ∆ ⇓ {ν 6 }; -after a new increasing sequence, we get (0 ≤ j ≤ i ≤ j + 3); -the decreasing sequence converges in one step, giving the best possible polyhedral invariant (0 ≤ j ≤ i ≤ j + 3, j ≤ 4). Notice that the constraint i ≤ j + 3 doesn't appear in the program, so it could not be chosen as a "threshold" for widening.

End of increasing sequence End of decreasing sequence Restarting at 0 ≤ i ≤ 3, j = 0

End of increasing sequence

End of decreasing sequence 5 Some improvements

0 ≤ j ≤ i ≤ j + 3 0 ≤ j ≤ i ≤ j + 3, j ≤ 4 0 ≤ j ≤ i 0 ≤ j ≤ i,
Of course, our method involves the computation of new iteration sequences, which may look expensive. We propose here two simple improvements to limit the cost of this computation. The first one is obvious: since the result, in general, is the greatest lower bound of the classical solution Z ∆ and the new limit V ∆ , one can intersect with Z ∆ each term of the new increasing sequence. In our Example 3, this would force the convergence of the new increasing sequence directly on the fixpoint (0 ≤ j ≤ i ≤ j + 3, j ≤ 4), without need of a decreasing sequence.

The second improvement is a compromise: it saves computation, but may lose precision (although we did not find any example where it happens). The selected components {Z ∆ j | ν j ∈ S} are intended to lead to improvements of other components, but in general, they are not supposed to be improved by the new iteration. With this idea in mind, we can limit the computation of the new iteration to the part of the CFG which doesn't influence only the components in S. More precisely, a component at node ν k must be computed again only if there is a path from ν k to a widening node which doesn't intersect S. In our Example 3, this would limit the new iteration to the subgraph {ν 7 , ν 9 , ν 11 }, since any path from {ν 8 , ν 10 , ν 3 , ν 4 } to the widening node ν 7 goes through the selected node ν 6 .

Experimental results

Our technique has been implemented inside our prototype static analyser, called Pagai, which computes numerical invariants in programs expressed in the LLVM internal representation [START_REF] Lattner | LLVM: a compilation framework fopr lifelong program analysis & transformation[END_REF]. In this representation, a function is a graph of basic blocks. The analyser takes as input such an LLVM file (that can be obtained from a C, C++, Fortran program by llvm-gcc or clang), and outputs for each basic block a numerical inductive invariant over the variables that are live at the head of this block. We can choose among several abstract domains : convex polyhedra, octagons, intervals, etc. through the Apron library [START_REF] Jeannet | Apron: a library of numerical abstract domains for static analysis[END_REF]. Since Pagai is an intra-procedural analyser, we can apply function inlining to obtain more precise results. We can also apply some LLVM optimisation passes, such as loop unrolling, or promoting memory variables to registers (mem2reg), in order to increase precision.

In the current state of our implementation, we only choose as improving components basic-blocks that are direct predecessors of a widening point, i.e., we don't apply our criteria C0 and C1 in their full generality. We apply only the first improvement proposed in §5: during the new increasing sequence, we intersect our new result with Z ∆ at each step.

We compared our technique with the classical abstract interpretation with standard widening/narrowing, on a variety of benchmarks.

The benchmark from the Mälardalen WCET research group1 contains interesting programs such as sorts, matrix transformations, fft, etc. These programs have been instrumented with a variable that counts the number of instructions being executed. These 98 functions have been analysed, using the polyhedra abstract domain. For 69 functions, the results of the new method are the same as the classical one, with a negligible time overhead (1.008 factor). For the other 29 functions, the new methods gives better results at 35% of widening points, with a 1.76 overhead factor. So, on this benchmark, not only the results are better for a significant subset of functions, but the new method costs almost nothing when it doesn't improve the results.

However, these encouraging experimental conclusions should not be overestimated: on a benchmark made of various highly used GNU functions (e.g., a2ps, gawk, gnuchess, gnugo, grep, gzip, lapack, make, sed and tar), the results are improved only at 4.14% of the widening points, with an overhead factor of 1.56 even on non improved functions.

A claim of the present paper is that the information about the end of executions can be as rich, complex and useful than the one derived from their beginning. To permit a better gathering of this information, we presented a method to improve the solution obtained by classical analysis: this solution is projected on some of its components, and the result of the projection is used to start a new pair of increasing and decreasing sequences.

The method is independent of the abstract lattice, it is compatible with any smart widening operator and any iteration strategy.

  A classical example where the decreasing sequence reaches the least fixpoint A nested loop prevents the decreasing sequence to get precise results

Fig. 3 .

 3 Fig. 3. CFG of the example 1.b

  Fig. 4. Example 3

Fig. 5 .

 5 Fig. 5. Analysis of Example 3

www.mrtc.mdh.se/projects/wcet/benchmarks.html

Acknowledgement: We are indebted to Laure Gonnord and David Monniaux for having put our attention on the problem -and specially on Examples 1 and 2 -and for helpful discussions.

⋆ This work has been partially supported by the ASOPT project of the "Agence Nationale de la Recherche" of the French Ministry of Research ⋆⋆ Verimag is a joint laboratory of Université Joseph Fourier, CNRS and Grenoble-INP.


V Monchiet
email: vincent.monchiet@univ-paris-est.fr

G Bonnet

Keywords: Sixth order tensor, Anisotropy, Inhomogeneity problem, Cylindrical inclusions

In this paper we provide a complete and irreducible representation for transversely isotropic sixth order tensors having minor symmetries. Such tensors appear in some practical problems of elasticity for which their inversion is required. For this kind of tensors, we provide an irreducible basis which possesses some remarkable properties, allowing us to provide a representation in a compact form which uses two scalars and three matrices of dimension 2, 3 and 4. It is shown that the calculation of sum, product and inverse of transversely isotropic sixth order tensors is greatly simplified by using this new formalism and appears to be appropriate for deriving new various solutions to some practical problems in mechanics which use such kinds of higher order tensors. For instance, we derive the fields within a cylindrical inhomogeneity submitted to remote gradient of strain. The method of resolution uses the Eshelby equivalent inclusion method extended to the case of a polynomial type eigenstrain. It is shown that the approach leads to a linear system involving a sixth order tensor whose closed form solution is derived by means of the tensorial formalism introduced in the first part of the paper.

Introduction

The characterization and classification of the material symmetries 1 for an anisotropic fourth order tensor have been investigated by several researchers, notably by Love [START_REF] Love | A Treatise on the Mathematical Theory of Elasticity, Fourth Edition[END_REF], Voigt [START_REF] Voigt | Lehrbuch der Krystallphysik[END_REF], Gurtin [START_REF] Gurtin | The Linear Theory of Elasticity[END_REF] and Thurston [START_REF] Thurston | Waves in solids[END_REF], Nye [START_REF] Nye | Physical properties of crystals[END_REF], Hou and Del Pierro [START_REF] Hou | On the completeness of the crystallographic symmetries in the description of the symmetries of the elastic tensor[END_REF], Forte and Vianello [START_REF] Forte | Symmetry classes for elasticity tensors[END_REF], Chadwick et al. [START_REF] Chadwick | A new proof that the number of linear elastic symmetries is eight[END_REF] among many others. In these studies, the authors express the different forms of the elasticity tensor into a matrix of dimension 6 × 6 (by means of Voigt type notations) in the cartesian system. In the paper of Walpole [START_REF] Walpole | Fourth-rank tensors of the thirthy-two crystal classes: multiplication table[END_REF] and later in [START_REF] Bona | Coordinate-free characterization of symmetry classes of elasticity tensors[END_REF], coordinate free representations of the elasticity tensors are derived by means of irreducible bases which prove to have interesting properties for the tensorial calculus. Indeed, by inspecting the properties of the multiplication table associated to these irreducible bases, symbolic representations which use scalars and matrices have been provided and shown to be appropriate for performing the standard tensorial operations such as the multiplication between two tensors and the inversion of a fourth order tensor. These "symbolic" or "compact" representations are very useful and suitable for deriving the solutions of elasticity problems for which the multiplication or more specially the inversion of fourth order tensors are required. Some fundamental theories as well as practical problems in the area of mechanics require the introduction of tensors whose order is higher than 4. For instance, the theories of gradient elasticity or the theory of microstructures in elasticity (see Toupin [START_REF] Toupin | Elastic materials with couple-stresses[END_REF], Mindlin [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], Mindlin & Eshel [START_REF] Mindlin | On first strain-gradient theory in linear elasticity[END_REF], Suiker and Chang [START_REF] Suiker | Application of higher order theory for formulating enhanced continuum models[END_REF]) introduce gradient elastic tensors of order 5, 6, ... in addition to the commonly used fourth order elasticity tensor. The use of such high order tensors raises a number of fundamental questions, especially concerning their inversion or the condition of positiveness when computing the elastic energy from these tensors. On the other hand, solutions to special elasticity problems, such as the solution to higher order inhomogeneity problems, require the inversion of tensors of order 6, 8, etc. (see Sendeckyj [START_REF] Sendeckyj | Ellipsoidal Inhomogeneity Problem[END_REF], Moschovidis [START_REF] Moschovidis | Two ellipsoidal inhomogeneities and related problems treated by the equivalent inclusion method[END_REF], Moschovidis and Mura [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], Asaro and Barnett [START_REF] Asaro | The non-uniform transformation strain problem for an anisotropic ellipsoidal inclusion[END_REF], Furuhashi and Mura, [START_REF] Furuhashi | On the equivalent inclusion method and impotent eigenstrains[END_REF], Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]). Although the tensorial equation produced in this case can be solved numerically, it is obviously of interest to derive closed form solutions which can cover many applications. Inverting such tensors involves the algebra of these higher order tensors. On the other hand, an irreducible basis for isotropic sixth and eighth order tensors having the minor symmetries 2 and a closed-form solution of higher order inhomogeneity problems have been recently provided by Monchiet and Bonnet [START_REF] Monchiet | Inversion of higher order isotropic tensors with minor symmetries and solution of higher order heterogeneity problems[END_REF]. However, the hypothesis of isotropy restricts the applicability of this new formalism for handling the problem of a spherical inhomogeneity. In this paper we propose to extend these results to the case of transversely isotropic sixth order tensors. We derive an irreducible basis, constituted of 31 elements, for any transversely isotropic sixth order tensors. The proposed representation is coordinate free since all these tensors are constructed as the outer products of elementary tensors attached to the direction of transverse isotropy. The multiplication tables are provided in the paper as well as a compact representation which uses two scalars and three matrices of dimension 2, 3 and 4. The efficiency of this new formalism and its ability to derive closed form solution to linear elasticity problems is thereafter illustrated in the case of a cylindrical inhomogeneity subjected to a remote uniform gradient of strain.

The paper is organized as follows:

-Section 2 states about the definition of transversely isotropic tensors and addresses the problem related to their inversion. -In section 3, we prove that the dimension of a transversely isotropic sixth order tensor is 31. -Section 4 then proposes an irreducible basis for such tensors and provides the multiplication tables as well as the compact representation. -Section 5 expounds the problem of an ellipsoidal inhomogeneity subjected to a gradient of strain and its method of resolution based on the Eshelby equivalent inclusion problem. Closed form solutions are then provided for the cylindrical inhomogeneity.

Below are provided some specific notations used in the paper:

a tensor of rank 1, a
tensor of rank 2 and 3, A tensor of rank 4 and 6, ⊗ tensorial product between two tensors, ⊙ n n th contraction between two tensors (ex: (A ⊙ 3 a) ijk = A ijkpqr a pqr ).

Statement of the problem

Consider the following linear relation:

A ⊙ n a = b (1)
involving two tensors a and b of order n and a 2n th order tensor A. In (1), the symbol "⊙ n " represents the n th contraction between two tensors such that (A ⊙ n a) i..j = A i..jp..q a p..q . In this paper we only examine the cases n = 2 and n = 3.

In the case n = 2, tensors a and b are symmetric second order tensors, a ij = a ji and b ij = b ji . Tensor A is a fourth order tensor having the minor symmetries,

A ijpq = A jipq = A ijqp , but not the major symmetry A ijpq = A pqij .
For the case n = 3, tensors a and b are third order tensors which are symmetric according to their first two indices: a ijk = a jik and b ijk = b jik . Due to the symmetries of a and b, A ijkpqr is invariant by any permutation of indices (i, j), and (p, q); these symmetries are called minor symmetries for the sixth order tensor. A tensor of components A i...j in the cartesian frame (e 1 , e 2 , e 3 ) is assumed symmetric with respect to the linear transformation R if:

A i...j = R ip ...R jq A p...q . (2)
We assume that: • A is invariant by the symmetry with respect to the plane Ox 1 x 2 characterized by the transformation of components R ij = δ ij -2n i n j and n is the vector e 3 .

• A i..j is invariant through any rotation by an angle of θ around the axis Ox 3 . In other words, A ijkpqr remains unchanged under the transformation [START_REF] Benedikt | On elastic interactions between spherical inclusions by the equivalent inclusion method[END_REF] where the components of the matrix R ij are given by:

R ij =      cos(θ) sin(θ) 0 -sin(θ) cos(θ) 0 0 0 1      . ( 3
)
With the above properties of invariance, tensor A is transversely isotropic and we denote by E 2n the space of such tensors (with n = 2 for fourth order tensors and n = 3 for sixth order tensors).

For the case n = 2, tensor A is defined by 6 independent coefficients. Providing a basis in the case of such symmetric tensors has been effected in [START_REF] Bona | Coordinate-free characterization of symmetry classes of elasticity tensors[END_REF]. However, the product of two such symmetric tensors is not generally symmetric and it is therefore important to study simultaneously symmetric and unsymmetric tensors. So, an irreducible basis for such tensors has been first provided in [START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF] and later in [START_REF] Walpole | Fourth-rank tensors of the thirthy-two crystal classes: multiplication table[END_REF]. In the following, we recall the principal results about this decomposition. Any fourth order tensor A ∈ E 4 can be constructed by the linear combination:

A = n=6 n=1 a n T n (4)
where the components of tensors T n for n = 1..6 are given by:

                             [T 1 ] ijpq = (π ip π jq + π iq π jp -π ij π pq )/2, [T 2 ] ijpq = (π ip n j n q + π iq n j n p + π jp n i n q + π jq n i n p )/2, [T 3 ] ijpq = π ij π pq /2, [T 4 ] ijpq = π ij n p n q / √ 2, [T 5 ] ijpq = n i n j π pq / √ 2, [T 6 ] ijpq = n i n j n p n q (5)
where π ij is defined by:

π ij = δ ij -n i n j (6)
and n i are the components of the normal unit vector n defining the direction orthogonal to the plane of transverse isotropy. The double contraction between two tensors taken from (T 1 , .., T 6 ) is provided in table 1.

⊙ 2 T 1 T 2 T 3 T 4 T 5 T 6 T 1 T 1 0 0 0 0 0 T 2 0 T 2 0 0 0 0 T 3 0 0 T 3 T 4 0 0 T 4 0 0 0 0 T 3 T 4 T 5 0 0 T 5 T 6 0 0 T 6 0 0 0 0 T 5 T 6
Table 1: The double contraction between the T n for n = 1..6

It can be observed that the set of tensors {T 1 , T 2 , T 3 , T 4 , T 5 , T 6 } is constituted of three groups which are independents for the composition ⊙ 2 , they are: {T 1 }, {T 2 }, {T 3 , T 4 , T 5 , T 6 }. Note that tensor A ∈ E 4 has the major symmetry only if a 4 = a 5 in (4). The table shows clearly that the product of two symmetric tensors is generally not symmetric. For completeness, we provide the relations giving the coefficients a 1 , .., a 6 in (4) as function of the components of tensor A:

a 1 = A ijpq π ip π jq /2 -A ijpq π ij π pq /4, a 2 = A ijpq π ip n j n q , a 3 = A ipjq π ip π jq /2, a 4 = A ijpq π pq n i n j / √ 2, a 5 = A ijpq π ij n p n q / √ 2, a 6 = A ijpq n i n j n p n q . ( 7
)
Denoting T n (a) = T n ⊙ 2 a for n = 1..6, their components are:

                             T 1 ij (a) = π ip π jq a pq -π ij a pq π pq /2, T 2 ij (a) = (π ip n j + π jp n i )a pq n q , T 3 ij (a) = π ij a pq π pq /2, T 4 ij (a) = π ij a pq n p n q / √ 2, T 5 ij (a) = n i n j a pq π pq / √ 2, T 6 ij (a) = n i n j a pq n p n q . ( 8
)
It follows that the linear relation [START_REF] Asaro | The non-uniform transformation strain problem for an anisotropic ellipsoidal inclusion[END_REF] for n = 2, can be rewritten into the form:

               T 1 (a) = a 1 T 1 (a), T 2 (a) = a 2 T 2 (a), T 3 (a) = a 3 T 3 (a) + a 4 T 4 (a), T 6 (a) = a 5 T 5 (a) + a 6 T 6 (a) (9)
where the a n for n = 1..6 are the components of tensor A along the basis (T 1 , ..., T 6 ). Owing to these results, it appears that the components of tensor A can be represented by two scalars and a 2 × 2 matrix. This compact representation is written into the symbolic form: A = {a 1 , a 2 , A} where the matrix A is defined by:

A =   a 3 a 4 a 5 a 6   . ( 10
)
It turns out that this symbolic notation A = {a 1 , a 2 , A} is useful for performing the standard tensorial operations. For instance, let us introduce a second tensor B, we denote by b 1 , ..b 6 its components along the basis (T 1 , ..., T 6 ) and {b 1 , b 2 , B} its compact representation. The product between A and B defines the fourth order tensor C = A⊙ 2 B whose components are

C = {c 1 = a 1 b 1 , c 2 = a 2 b 2 , C = A.
B} where A.B represents the standard matrix product between A and B. The inverse of tensor A is simply obtained by inverting the elements of its compact representation: A -1 = {1/a 1 , 1/a 2 , A -1 } where A -1 denotes the inverse of the 2 × 2 matrix. Such a representation is useful for doing the classical tensorial operations and is convenient for many practical problems in mechanics such as Eshelby's inhomogeneity problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF]. In the following section, we aim at deriving a similar representation for transversely isotropic sixth order tensors.

3 Dimension of the space containing transversely isotropic and symmetric sixth order tensors

A key point in determining the representation of a tensor of a specific type by an irreducible basis lies in the calculation of the total number of elements of the basis. In this section, we address a general method for computing the dimension of the space containing sixth order tensors having the minor symmetries and we provide results for the particular case of the transversely isotropic symmetry.

As effected for the case of fourth order tensors, Voigt type notations are suitable for representing such tensors by a 6×6 matrice. In the following, we extend it to the case of a sixth order tensor. This is achieved by a vector representation of third-order tensors. Any third order tensor a, which is symmetric according to its two first indices (a ijk = a jik ) depends on 18 independent components and can the been represented along tensors v n for n = 1..18 which are defined by :

                             v 1 = f 1 ⊗ e 1 , v 2 = f 2 ⊗ e 1 , v 3 = f 3 ⊗ e 1 , v 4 = f 6 ⊗ e 2 , v 5 = f 5 ⊗ e 3 , v 6 = f 2 ⊗ e 2 , v 7 = f 1 ⊗ e 2 , v 8 = f 3 ⊗ e 2 , v 9 = f 6 ⊗ e 1 , v 10 = f 4 ⊗ e 3 , v 11 = f 3 ⊗ e 3 , v 12 = f 1 ⊗ e 3 , v 13 = f 2 ⊗ e 3 , v 14 = f 5 ⊗ e 1 , v 15 = f 4 ⊗ e 2 , v 16 = f 4 ⊗ e 1 , v 17 = f 5 ⊗ e 2 , v 18 = f 6 ⊗ e 3 (11)
where e 1 , e 2 , e 3 are the unit vectors of the cartesian basis while the second order tensors f n for n = 1..6 are defined by :

     f 1 = e 1 ⊗ e 1 , f 2 = e 2 ⊗ e 2 , f 3 = e 3 ⊗ e 3 , f 4 = √ 2e 2 ⊗ s e 3 , v 5 = √ 2e 1 ⊗ s e 3 , f 6 = √ 2e 1 ⊗ s e 2 (12)
where ⊗ s represents the symmetrized tensorial product of two vectors:

u⊗ s v = (u ⊗ v + v ⊗ u)/2.
With the above definitions for tensors v n , the basis v 1 ..v 18 is orthonormal for the composition "⊙ 3 " since :

v r ⊙ 3 v s = δ rs (13)
for any r, s = 1..18. In relation (1) (for the case n = 3), tensors a, b and A are decomposed along the basis v n for n = 1..18 :

a = a I v I , b = b I v I , A = A IJ v I ⊗ v J . (14)
By convention the upper case indices vary from 1 to 18. Einstein summation convention for repeated indices is applied in the above expressions. For instance a I v I = a 1 v 1 + ... + a 18 v 18 . The third order tensors a and b can then be replaced by two vectors of dimension 18 having respectively the components a I and b I for I = 1..18. In the linear relation b = A ⊙ 3 a, the components of b depend on the ones of a by: b I = A IJ a J where A IJ is a matrix of dimension 18 × 18 which depends on 324 independent coefficients. When tensor A has the major symmetry, matrix A IJ is symmetric and depends on 171 independent coefficients. Tensors v I ⊗ v J for I, J = 1..18 constitute then an irreducible basis for general sixth order tensors. When A is invariant by rotation or reflection (cf. relation ( 2)), the total number of independent coefficients of A can be strongly reduced and the definition of a new irreducible basis is required.

By adopting the above contracted notations, the transformation rule (2) reads:

A IJ = R IP R JQ A P Q (15)
where matrix R IJ are the components of tensor R defined by R ijkpqr = (R ip R jq + R iq R jp )R kr /2 with the above contracted notations. In the case of a plane symmetry with respect to the plane Ox 1 x 2 , R IJ is diagonal and is given by:

R IJ = diag(1, 1, .., 1 10 times , -1, -1, .., -1 8 times
). ( 16)

In the case of a rotation of an angle θ around the axis Ox 3 , the matrix R IJ is obtained by the concatenation of different matrices as follows:

R IJ =      R 1 ij R 2 ij 0 -R 2 ij R 1 ij 0 0 0 R 3 ij      . (17)
In the above expressions, matrices R 1 ij , R 2 ij are squared matrices of dimension 5, defined by:

R 1 ij = p          p 2 q 2 0 √ 2q 2 0 q 2 p 2 0 - √ 2q 2 0 0 0 1 0 0 √ 2q 2 - √ 2q 2 0 p 2 -q 2 0 0 0 0 0 1          , R 2 ij = q           q 2 p 2 0 √ 2p 2 0 p 2 q 2 0 - √ 2p 2 0 0 0 1 0 0 √ 2p 2 - √ 2p 2 0 q 2 -p 2 0 0 0 0 0 1           . (18)
R 3 ij is squared matrix of dimension 8; its components are:

R 3 ij =                   1 0 0 0 0 0 0 0 0 p 2 q 2 0 0 0 0 - √ 2pq 0 q 2 p 2 0 0 0 0 √ 2pq 0 0 0 p 2 q 2 -pq -pq 0 0 0 0 q 2 p 2 pq pq 0 0 0 0 -pq pq p 2 -q 2 0 0 0 0 -pq pq -q 2 p 2 0 0 √ 2pq - √ 2pq 0 0 0 0 p 2 -q 2                   . ( 19
)
The invariance condition ( 15) is now verified with expressions ( 16) and ( 17) for the matrix R IJ . These computations are elementary and lead to the following form for A IJ :

A IJ =         f ij 0 0 0 0 f ij 0 0 0 0 g ij 0 0 0 0 h ij         (20)
where f ij , g ij and h ij are squared matrices of dimension 5, 5 and 3 respectively. The following equalities hold for the components of these matrices:

                                   f 21 -f 12 = √ 2(f 14 + f 41 ), f 11 -f 22 = √ 2(f 24 + f 41 ) = √ 2(f 42 + f 14 ), 2f 44 = f 11 + f 22 -f 21 -f 12 , √ 2f 45 = f 15 -f 25 , √ 2f 43 = f 13 -f 23 , √ 2f 54 = f 51 -f 52 , f 34 = f 31 -f 32 ,
g 12 = g 13 , g 14 = g 15 , g 21 = g 31 , g 22 = g 33 , g 23 = g 32 , g 24 = g 35 , g 25 = g 34 , g 41 = g 51 , g 42 = g 53 , g 43 = g 52 , g 44 = g 55 , g 45 = g 54 ,

h 22 = h 11 , h 23 = h 13 = g 42 -g 43 , √ 2h 32 = √ 2h 31 = g 24 -g 34 , h 21 = h 12 , h 33 = g 22 -g 23 , h 11 + h 12 = g 44 -g 45 . (21)
Matrices f ij and g ij for i, j = 1..5 and h ij for i, j = 1..3 are defined by 59 coefficients. Moreover, by considering the above 28 equalities, it follows that the total number of independent coefficients for matrix A IJ is 31.

An irreducible basis for transversely isotropic sixth order tensors

In this section we propose a canonical basis for transversely isotropic sixth order tensors. The dimension of A ∈ E 6 being 31, an irreducible basis of such kind of tensor comprises 31 independent tensors. The method of construction uses a decomposition of a third order tensor related to the transverse isotropy.

First, we decompose the third order tensor a, of components a ijk , on the form:

a ijk = w ijk + u ij n k + n j v ik + n i v jk + n i n j r k +(n i s j + n j s i )n k + γn i n j n k (22)
where n i are the components of the vector n of transverse isotropy. The inner product between w ijk , u ij , v ij , r i , s i with n i is null:

w ijk n k = w ikj n k = w kij n k = 0, u ij n j = u ji n j = v ij n j = v ij n i = 0, r i n i = s i n i = 0. ( 23
)
The computation of w ijk , u ij , v ij , r i , s i and γ from equations ( 22) and ( 23) leads to the following expressions:

                             γ = a ijk n i n j n k , r i = π ir a pqr n p n q , s i = π ip a pqr n q n r , u ij = π ip π jq a pqr n r , v ij = π iq π jr a pqr n p , w ijk = π ip π jq π kr a pqr ( 24
)
where it is recalled that

π ij = δ ij -n i n j . Note that u ij is symmetric, u ij = u ji ; v ij which is non-symmetric,
is decomposed into a symmetric and a skew symmetric tensor:

v ij = v s ij + v a ij (25)
with:

v s ij = 1 2 (v ij + v ji ), v a ij = 1 2 (v ij -v ji ). ( 26
)
Tensors u and v s are decomposed into two-dimensional spherical and deviatoric parts:

u ij = u ij -π ij α, v s ij = v s ij -π ij β (27)
with:

α = 1 2 π pq a pqr n r , β = 1 2 π qr a pqr n p . (28)
Finally we introduce a similar decomposition for w ijk :

w ijk = w ijk + p i π jk + p j π ik + q k π ij (29)
with:

p i = w ijk π jk /2 -w jki π jk /4, q i = 3w jki π jk /4 -w ijk π jk /2 (30)
where the third order tensor w is traceless: w ijj = w jji = 0 and is invariant by any permutation of its indices, w ijk = w kij = w jki = w kji = w ikj = w jik . Let us recall that any third order tensor having the minor symmetries depends on 18 independent coefficients. Each one of vectors p, q, r and s, of components p i , q i , r i and s i , depends on 2 coefficients. The second order tensors u, v s and v a are defined by 2, 2 and 1 independent coefficient. The third order tensor w depends on 2 independent coefficients. The total number of independent coefficients of each element of the decomposition are recalled below:

α β γ p q r s u v s v a w a 1 1 1 2 2 2 2 2 2 1 2 18
Table 1 Number of independent coefficients of the elements of the decomposition of the third order tensor a

The various outer products of w ijk , v s ij , u ij , v a ij , t i , s i , x, v m and u m with π ij and n i provide 31 third order tensors having the following components:

T ijk (a) = w ijk , T 2 ijk (a) = n i v a jk + n j v a ik , T ijk (a) = u ij n k , T 4 ijk (a) = 2v s ij n k , T ijk (a) = (n i u jk + n j u ik )/2, T 6 ijk (a) = n i v s jk + n j v s ik , T ijk (a) = απ ij n k , T 8 ijk (a) = 2βπ ij n k , T ijk (a) = γπ ij n k , T 10 ijk (a) = α(n i π jk + n j π ik )/2, T ijk (a) = β(n i π jk + n j π ik ), T 12 ijk (a) = γ(n i π jk + n j π ik )/2, T ijk (a) = αn i n j n k , T 14 ijk (a) = 2βn i n j n k , T ijk (a) = γn i n j n k , T 16 ijk (a) = p i π jk + p j π ik , T ijk (a) = (q i π jk + q j π ik )/2, T 18 ijk (a) = (r i π jk + r j π ik )/2, T ijk (a) = (s i π jk + s j π ik )/2, T 20 ijk (a) = 2π ij p k , T ijk (a) = π ij q k , T 22 ijk (a) = π ij r k , T ijk (a) = π ij s k , T 24 ijk (a) = 2n i n j p k , T ijk (a) = n i n j q k , T 26 ijk (a) = n i n j r k , T ijk (a) = n i n j s k , T 28 ijk (a) = 2(p i n j + p j n i )n k , T ijk (a) = (q i n j + q j n i )n k , T 30 ijk (a) = (r i n j + r j n i )n k , T ijk (a) = (s i n j + s j n i )n k . (31)
In relation b = A ⊙ 3 a, the third tensor b is a transversely isotropic and linear function of the third order tensor a. The third order tensors T n (a) for n = 1..31 provide a complete and irreducible representation for the third order valued tensor function b. This is proved in appendix A. We now introduce the sixth order tensors T n for n = 1..31 such that T n (a) = T n ⊙ 3 a. The components of these tensors are given below:

[T 1 ] ijkpqr = π kr (π ip π jq + π jp π iq )/2 -3π ij π pq π kr /4 +π ij (π pr π kq + π qr π kp )/4 + π pq (π ir π jk + π jr π ik )/4 -π ik (π jp π qr + π jq π pr )/4 -π jk (π ip π qr + π iq π pr )/4, [T 2 ] ijkpqr = (n i n p π jq + n i n q π jp + n j n p π iq + n j n q π ip )π kr /4 -(π ir n j + π jr n i )n q π kp /4 -(π ir n j + π jr n i )n p π kq /4, [T 3 ] ijkpqr = (π ip π jq + π iq π jp -π ij π pq )n k n r /2, [T 4 ] ijkpqr = -(n p π qr + n q π pr )n k π ij /2 +(π iq π jr + π jq π ir )n p n k /2 + (π ip π jr + π jp π ir )n q n k /2, [T 5 ] ijkpqr = -(n i π jk + n j π ik )π pq n r /4 +(π ip π kq + π iq π kp )n j n r /4 + (π jp π kq + π jq π kp )n i n r /4, [T 6 ] ijkpqr = (n i n p π jq + n i n q π jp + n j n p π iq + n j n q π ip )π kr /4 +(π ir n j n q + π jr n i n q )π kp /4 + (π ir n j n p + π jr n i n p )π kq /4 -(n i n p π qr + n i n q π pr )π jk /4 -(n j n p π qr + n j n q π pr )π ik /4, [T 7 ] ijkpqr = π ij π pq n k n r /2, [T 8 ] ijkpqr = (n p π qr + n q π pr )n k π ij /2, [T 9 ] ijkpqr = π ij n p n q n k n r , [T 10 ] ijkpqr = (n i π jk + n j π ik )n r π pq /4, [T 11 ] ijkpqr = (n p π qr + n q π pr )(n i π jk + n j π ik )/4, [T 12 ] ijkpqr = (π ik n j + π jk n i )n r n p n q /2, [T 13 ] ijkpqr = n i n j π pq n k n r /2, [T 14 ] ijkpqr = (π pr n q + π qr n p )n i n j n k /2, [T 15 ] ijkpqr = n i n j n k n p n q n r , [T 16 ] ijkpqr = π jk (π ip π qr + π iq π pr -π pq π ir )/4 +π ik (π jp π qr + π jq π pr -π pq π jr )/4, (32)
[T 17 ] ijkpqr = π jk (3π pq π ir -π ip π qr -π iq π pr )/8 +π ik (3π pq π jr -π jp π qr -π jq π pr )/8, [T 18 ] ijkpqr = (π ir π jk + π jr π ik )n p n q /2,

[T 19 ] ijkpqr = (π ip n q + π iq n p )π jk n r /4 + (π jp n q + π jq n p )π ik n r /4, [T 20 ] ijkpqr = π ij (π pr π kq + π qr π kp -π pq π kr )/2, [T 21 ] ijkpqr = π ij (3π pq π kr -π pr π kq -π qr π kp )/4, [T 22 ] ijkpqr = π ij n p n q π kr , [T 23 ] ijkpqr = (n p π qk + n q π pk )n r π ij /2, [T 24 ] ijkpqr = (π pr π qk + π qr π pk -π pq π kr )n i n j /2, [T 25 ] ijkpqr = (3π pq π kr -π pr π qk -π qr π pk )n i n j /4, [T 26 ] ijkpqr = n i n j n p n q π kr , [T 27 ] ijkpqr = (π pk n q + π qk n p )n i n j n r /2, [T 28 ] ijkpqr = (π ip π qr + π iq π pr )n j n k /2 + (π jp π qr + π jq π pr )n i n k /2 -(n i π jr + n j π ir )n k π pq /2, [T 29 ] ijkpqr = 3(n i π jr + n j π ir )n k π pq /4 -(π ip π qr + π iq π pr )n j n k /4 -(π jp π qr + π jq π pr )n i n k /4, [T 30 ] ijkpqr = (π ir n j + π jr n i )n k n p n q , [T 31 ] ijkpqr = (π ip n j n q + π iq n j n p + π jp n i n q + π jq n i n p )n k n r /2. ( 33
)
The set of tensors T n for n = 1...31 constitutes then an irreducible basis for any transversely isotropic sixth order tensor having the minor symmetries.

Then any tensor A ∈ E 6 can be read:

A = n=31 n=1 a n T n (34)
where the a n are the components of tensor A in the basis (T It can be noticed that the irreducible basis for fourth order tensors, defined in [START_REF] Forte | Symmetry classes for elasticity tensors[END_REF], is also constituted of idempotent tensors (T 1 , T 2 , T 3 and T 6 ) and nilpotent tensors (they are T 4 and T 5 ). Note also that similar structures has been also obtained for sixth and eighth order isotropic tensors in Monchiet and Bonnet [START_REF] Monchiet | Inversion of higher order isotropic tensors with minor symmetries and solution of higher order heterogeneity problems[END_REF], [START_REF] Monchiet | On the inversion of non symmetric sixth order isotropic tensors and conditions of positiveness of third order tensor valued quadratic functions[END_REF]. It can be observed that the set (E 6 , ⊙ 3 , I) defines a monoid (an algebraic structure with a single associative binary operation and an identity element). I represents the identity for isotropic sixth order tensors having the minor symmetries and is such that I ⊙ 3 a = a for any third order tensor a. The components of the sixth order identity tensor are I ijkpqr = I ijpq δ kr where I ijpq = (δ ip δ jq + δ iq δ jp )/2 are the components of the identity tensor for symmetric fourth order tensors. The decomposition of I along the basis (T

⊙ 3 T 3 T 4 T 5 T 6 T 3 T 3 T 4 0 0 T 4 0 0 T 3 T 4 T 5 T 5 T 6 0 0 T 6 0 0 T 5 T 6
Table 3: The triple contraction between the T n for n = 3..6 where coefficients a n for n = 1..31 are the components of A as defined in (34). Consequently, it is possible to represent A by two scalars, and three matrices of dimension 2, 3 and 4 with the following symbolic notations:

                                                             T 1 (b) = a 1 T 1 (a), T 2 (b) = a 2 T 2 (a),
A = {a 1 , a 2 , A 1 , A 2 , A 3 }. ( 37
)
The components of A 1 , A 2 and A 3 are: and the a n for n = 1..31 are the components of A in the basis T n with n = 1..31 as defined in relation (34). For instance the identity I reads:

A 1 =   a 3 a 4 a 5 a 6   , A 2 =      a 7 a
I =                1, 1,   1 0 0 1   ,      1 0 0 0 1 0 0 0 1      ,         1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1                        . ( 39
)
Let us now introduce a second tensor B, whose components in the basis (T 1 , ..., T 31 ) are denoted by b n . We use the same compact representation for this tensor:

B = {b 1 , b 2 , B 1 , B 2 ,
c 1 = a 1 b 1 , c 2 = a 2 b 2 , C 1 = A 1 B 1 , C 2 = A 2 B 2 , C 3 = A 3 B 3 ( 40
)
where the standard matrix product rule is applied between the matrices A n and B n for obtaining the matrix C n . The inversion of tensor A is then obviously obtained by inverting all the elements of its compact representation:

A -1 = 1 a 1 , 1 a 2 , A -1 1 , A -1 2 , A -1 3 . ( 41
)
The existence of such an inverse is obviously dependent on the inverse of its components. As an illustration purpose, we now aim at applying this formalism for deriving closed-form solutions of the second order inhomogeneity problem.

Application to the second order inhomogenity problem

Eshelby's solutions [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF] for inclusions and for equivalent inhomogeneity problems are fundamental to many problems in material science, mechanics of composite, etc. In the terminology of Eshelby [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF] and Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF], an inclusion denotes a subdomain subjected to an eigenstrain while an inhomogeneity is a domain whose elastic properties differ from those of the surrounding medium. Eshelby's results [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF] are well known for the case of a prescribed constant eigenstrain in an infinite isotropic medium: it gives a constant strain field inside an ellipsoidal inclusion while the exterior point solution (outside the inclusion) is heterogeneous. The strain field within the inclusion is then explicitly given as a function of the eigenstrain by means of the well-known fourth order Eshelby's tensor. The Eshelby equivalent method handles the problem of a single ellipsoidal inhomogeneity by replacing it with an inclusion having properly chosen eigenstrains. The strain field within the inhomogeneity remains constant and can be expressed in terms of the remote constant strain field by inverting a fourth order tensor which does not possess the major symmetry. Later, Sendeckyj [START_REF] Sendeckyj | Ellipsoidal Inhomogeneity Problem[END_REF], Moschovidis [START_REF] Moschovidis | Two ellipsoidal inhomogeneities and related problems treated by the equivalent inclusion method[END_REF], Moschovidis and Mura [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], Asaro and Barnett [START_REF] Asaro | The non-uniform transformation strain problem for an anisotropic ellipsoidal inclusion[END_REF], generalize Eshelby's solution to the case of prescribed polynomial fields. In these studies, the following result has been proved: the strain in an ellipsoidal subdomain of an infinite linear elastic medium which undergoes an eigenstrain on the form of a polynomial of degree N , is also a polynomial with the same degree N . The expansions of the eigenstrain and of the interior point solution for the strain field along polynomial functions introduce tensors of order 2, 3, 4 etc and higher order Eshelby tensors of order 6, 8 etc. It has been pointed out in Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] that the strain disturbance due to a polynomial type remote strain field of degree N can be simulated by an appropriate polynomial eigenstrain field of degree N . Note that the utility of the result for N = 1 has been addressed by Eshelby [START_REF] Eshelby | Progress in solid mechanics[END_REF] for treating ellipsoidal inhomogeneities submitted to far-field torsion and flexure. Solutions to higher order inclusion problems are also the bases of numerous works dealing with multiple ellipsoidal inhomogeneities, among which Moschovidis and Mura [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF], Shodja and Sarvestani [START_REF] Shodja | Elastic fields in double inhomogeneity by the equivalent inclusion method[END_REF], Shodja et al. [START_REF] Shodja | Interacting cracks and ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], Benedikt et al. [START_REF] Benedikt | On elastic interactions between spherical inclusions by the equivalent inclusion method[END_REF].

In the next section we propose to recall the solution of the second order inclusion problem and we deal with the problem of an ellipsoidal inhomogeneity submitted to a remote gradient of strain. Closed form solutions are provided in section 5.2 for the case of a cylindrical inhomogeneity for which the sixth order tensor, which has to be inverted, presents a transverse isotropy, the axis of isotropy being parallel to the axis of the cylinder.

The general case of an ellipsoidal inhomogeneity

Consider an ellipsoidal inhomogeneity, of stiffness tensor C, embedded in an infinite elastic matrix of stiffness tensor C 0 . We denote by (x 1 , x 2 , x 3 ) the cartesian coordinates, (e 1 , e 2 , e 3 ) the orthonormal basis and by a 1 , a 2 , a 3 the radii of the ellipsoid along the three axes of the cartesian frame. The inclusion is located at the origin and its volume is defined by:

x 2 1 a 2 1 + x 2 2 a 2 2 + x 2 3 a 2 3 ≤ 1. (42)
The ellipsoidal inhomogeneity contains an isotropic elastic material, µ and λ being its Lamé coefficients. The infinite medium is also isotropic and its elastic moduli are denoted by µ 0 and λ 0 . This inhomogeneity is subjected to a remote strain field on the form ε ∞ ij = b ijk x k where x, of components x k , denotes the vector position. This problem is called "second order inhomogeneity problem" since an uniform gradient of strain is applied at the infinity instead of the constant strain considered by Eshelby [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF]. The strain disturbance due to the application of b ijk at infinity can be recovered by considering an appropriate inclusion problem. This auxiliary problem is the following: an ellipsoidal domain, defined by (42), is subjected to a prescribed eigenstrain on the form ε * ij (x) = e ijk x k (with e ijk = e jik ). As shown in Sendeckyj [START_REF] Sendeckyj | Ellipsoidal Inhomogeneity Problem[END_REF], Moschovidis [START_REF] Moschovidis | Two ellipsoidal inhomogeneities and related problems treated by the equivalent inclusion method[END_REF], the solution of this auxiliary problem is:

u i (x) = 1 8π(1 -ν 0 ) Ψ r,pqi e pqr -2ν 0 Φ r,i e ppr -4(1 -ν 0 )Φ r,k e ikr . (43)
The associated strain field, ε(x), reads:

ε ij (x) = 1 8π(1 -ν 0 ) Ψ r,pqij e pqr -2ν 0 Φ r,ij e ppr -2(1 -ν 0 )(Φ r,jk e ikr + Φ r,ik e jkr ) (44)
in which Φ i and Ψ i are respectively harmonic and bi-harmonic potentials whose expressions can be also found in the book of Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] (see equations 12.13 and 12.14) and are recalled in the appendix C. The strain field (44) for an interior point (i.e. within the ellipsoidal domain) is linear according to the vector position and can be put into the form:

ε ij (x) = S ijkpqr e pqr x k (45)
where S ijkpqr are the components of the sixth order Eshelby tensor denoted S. The associated stress field reads:

σ ij (x) = C 0 ijmn (S mnkpqr e pqr -e mnk )x k . ( 46
)
The components of the sixth order Eshelby tensor, S, have not been explicitly derived in Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]; their computation is detailed in appendix C and the final expression is given below:

8π(1 -ν 0 )S ijkpqr = δ ij δ kr δ pq T IP R + 2ν 0 I IR +2I ijpq δ kr T IJR + (1 -ν 0 )(I IR + I JR ) +2I pqkr δ ij T IKR + 2I ijkr δ pq T IJP + 2ν 0 I IJ +2(I ijpr δ kq + I ijqr δ kp ) T P QR + (1 -ν 0 )I KR +2(I ijkp δ qr + I ijkq δ pr ) T KP Q + (1 -ν 0 )I KR a 2 R . (47)
The following summation convention has been used: repeated lower case indices are summed from 1 to 3; upper case indices (unlike in section 3) take on the same values as the corresponding lower case ones but are not summed. For example, in the monomial a i a i b I , the repeated index is i and the upper case index takes the same value as i; it gives:

a i a i b I = a 2 1 b 1 +a 2 2 b 2 +a 2 3 b 3 .
The reader can refer to Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] for more details and examples about this summation convention. The I ij and T ijk coefficients are defined by the following elliptic integrals:

I ij = 2πa 1 a 2 a 3 +∞ 0 ds (a 2 i + s)(a 2 j + s)∆(s)
,

T ijk = 2πa 1 a 2 a 3 +∞ 0 sds (a 2 i + s)(a 2 j + s)(a 2 k + s)∆(s) (48)
with:

∆(s) = (a 2 1 + s) 1/2 (a 2 2 + s) 1/2 (a 2 3 + s) 1/2
. Note that I ij and T ijk are invariant by any permutation of their indices. The closed form expressions of the integrals appearing in (48) can be found in Moschovidis [START_REF] Moschovidis | Two ellipsoidal inhomogeneities and related problems treated by the equivalent inclusion method[END_REF], Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] for the case of a spheroidal inhomogeneity and particularly for the special case of a cylinder.

Let us now come back to the inhomogeneity problem. Using the equivalent Eshelby inclusion method, the eigenstrain ε * ij (x) is chosen as:

ε * ij (x) = -Y ijpq ε pq (x) (49)
where Y ijpq is defined by:

Y ijpq = 1 - k k 0 J ijpq + 1 - µ µ 0 K ijpq (50)
where k = λ + 2µ / 3 is the bulk modulus of the inclusion while k 0 = λ 0 + 2µ 0 /3 is the bulk modulus of the infinite matrix. In (50), J ijpq = δ ij δ pq /3 and K ijpq = I ijpq -J ijpq where it is recalled that I ijpq = (δ ip δ jq + δ iq δ jp )/2. The strain field in the ellipsoidal inhomogeneity is also linear according to the vector position and can then be written as ε ij (x) = a ijk x k where a ijk are the components of a constant third order tensor a which possesses the symmetry with respect to its two first indices. The strain field, solution of the inhomogeneity problem, is the sum of: (i) the prescribed remote strain field ε ∞ ij = b ijk x k , (ii) the disturbed strain field due to the applied eigenstrain ε * ij (x) = e ijk where e ijk is related to a ijk by e ijk = S ijkpqr a pqr . It follows that a is solution of a linear system on the form [START_REF] Asaro | The non-uniform transformation strain problem for an anisotropic ellipsoidal inclusion[END_REF] in which the components of A are given by:

A ijkpqr = I ijkpqr -S ijkmnr Y mnpq ( 51
)
where it is recalled that I ijkpqr are the components of the sixth order identity tensor having the minor symmetries (see section 4).

Closed-form solution for a cylindrical inhomogeneity

We propose here to consider the case of the cylindrical inhomogeneity. We put a 1 = a 2 = a (a being the radius of the circular section), while a 3 tends to infinity. The second order Eshelby tensor S defined in equation ( 47) with the definitions (48) can be represented along the basis (T 1 , ..., T 31 ) where e 3 is the direction of transverse isotropy. For the computation of the elements of the compact representation of S, denoted {s 1 , s 2 , S 1 , S 2 , S 3 }, we use the base change relations provided in appendix B. The final expressions for the components of the Eshelby tensor are (see appendix D for the detailed expressions of coefficients s n in the case of a spheroidal inhomogeneity):

s 1 = 1 4 3 -4ν 0 1 -ν 0 , s 2 = 0, S 1 = 1 4(1 -ν 0 )   3 -4ν 0 0 2(1 -2ν 0 ) 2(1 -ν 0 )   , S 2 = 1 2(1 -ν 0 )      1 0 2ν 0 2ν 0 2(1 -ν 0 ) 2 0 0 0      , S 3 = 1 4(1 -ν 0 )         4(1 -ν 0 ) 2 2ν 0 2(1 -2ν 0 ) 0 1 ν 0 1 0 0 0 4(1 -ν 0 ) 0 0 0 2(1 -ν 0 )         . ( 52
)
Tensor A, defined by (51), is also a transversely isotropic sixth order tensor. Its representation by means of the symbolic notation A = {a 1 , a 2 , A 1 , A 2 , A 3 } is then possible. For the computation of A, it is convenient to transform the inner product S ijkmnr Y mnpq into a triple contraction between two sixth order tensors: tensor S and tensor Y * whose components are Y * ijkpqr = Y ijpq δ kr where Y ijpq is defined by (50). In this manner, tensor A can be rewritten as A = I -S ⊙ 3 Y * where it is recalled that I is the identity for sixth order tensors having minor symmetries and for the triple contraction ⊙ 3 . The compact representation of tensor Y * is {ξ, ξ, Y 1 , Y 2 , Y 3 } with:

Y 1 =   ξ 0 0 ξ   , Y 2 =      (2κ + ξ)/3 0 (κ -ξ)/3 0 ξ 0 2(κ -ξ)/3 0 (κ + 2ξ)/3      , Y 3 =         ξ 0 0 0 (κ -ξ)/3 (2κ + ξ)/3 (κ -ξ)/3 0 (κ -ξ)/3 2(κ -ξ)/3 (κ + 2ξ)/3 0 0 0 0 ξ         (53)
in which κ and ξ are given by:

κ = 1 - k k 0 , ξ = 1 - µ µ 0 . (54)
Finally, the inversion of A is performed by inverting the elements of its compact representation. It leads to:

b 1 = 4(1 -ν 0 ) ∆ 1 , b 2 = 1, b 3 = b 1 , b 4 = 0, b 5 = 4ξ(1 -2ν 0 ) (2 -ξ)∆ 1 , b 6 = 2ξ 2 -ξ , b 7 = 6(1 -ν 0 ) ∆ 2 , b 8 = 0, b 9 = κ(1 + ν 0 ) -ξ(1 -2ν 0 ) ∆ 2 , b 10 = 2 2κ(1 + ν 0 ) -ξ(2 -ν 0 ) (1 -ξ)∆ 2 , b 11 = 1 1 -ξ , b 12 = 2κ(1 + ν 0 ) + 2ξ(2 -ν 0 ) -3κξ(1 + ν 0 ) (1 -ξ)∆ 2 , b 13 = b 14 = 0, b 15 = 1, b 16 = 12(1 -ν 0 ) -2κ(1 + ν 0 ) -ξ(1 -2ν 0 ) ∆ 3 , b 17 = 2 2κ(1 + ν 0 ) + ξ(1 -2ν 0 ) ∆ 3 , b 18 = 2κ(1 + ν 0 ) -2ξ(1 -2ν 0 ) ∆ 3 , b 19 = 6ξ(1 -2ν 0 ) ∆ 3 , b 20 = (κ -ξ)(1 + ν 0 ) ∆ 3 , b 21 = 2 6(1 -ν 0 ) -κ(1 + ν 0 ) + ξ(7ν 0 -5) ∆ 3 , b 22 = (κ(1 + ν 0 ) -ξ(1 -2ν 0 ))(1 -ξ) ∆ 3 , b 23 = 3ξ(1 -ξ)(1 -2ν 0 ) ∆ 3 - ξ 2 -ξ , b 24 = b 25 = 0, b 26 = 1, b 27 = 2ξ 2 -ξ , b 28 = b 29 = b 30 = 0, b 31 = 2 2 -ξ (55)
with:

∆ 1 = 4(1 -ν 0 ) + ξ(4ν 0 -3), ∆ 2 = 6(1 -ν 0 ) -2κ(1 + ν 0 ) -ξ(1 -2ν 0 ), ∆ 3 = 12(1 -ν 0 ) -4κ(1 + ν 0 ) + ξ(16ν 0 -11) + 2κξ(1 + ν 0 ) + ξ 2 (1 -2ν 0 ) (56)
which allows the explicit computation of the strain field inside the cylinder produced by the heterogeneity problem.

Conclusion

In this paper we have provided two new results:

• The first concerns the linear algebra of transversely isotropic sixth order tensors having the minor symmetries. For this class of symmetry, it has been shown that any sixth order tensor depends on 31 independent coefficients. An irreducible basis for such tensors has then been provided as well as a representation by means of two scalars and three matrices of dimensions 2, 3, and 4. This formalism is shown to be very useful for performing the classical tensorial operations since it reduces to matrix operations on matrices possessing small dimensions. • The second result has concerned the area of micromechanics of heterogenous materials. More precisely, we have derived the strain field occurring within a cylindrical elastic inhomogeneity embedded in an infinite elastic matrix and submitted to a remote gradient of strain. The Eshelby equivalent method, which has been employed for deriving the closed-form solution, requires the inversion of a transversely isotropic sixth order tensor for which the tensorial formalism introduced in the first part of the paper has been used.

The compact representation proposed in this paper for transversely isotropic sixth order tensor could be also used for deriving various results in the context of micromechanics of heterogeneous material (for instance, the problems related to the interaction between two inhomogeneities [START_REF] Shodja | Interacting cracks and ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF]). Moreover, this tensorial formalism can be also used for deriving the condition of positiveness of elastic energy of generalized continuum theories (as already done in the case of isotropic materials in [START_REF] Monchiet | On the inversion of non symmetric sixth order isotropic tensors and conditions of positiveness of third order tensor valued quadratic functions[END_REF]).

B Base change relations

Here we provide the expressions of coefficients an for n = 1..31, as defined in (34), as functions of the components A ijkpqr :

a 1 = A ijkpqr π ip π jq π kr /3 -

D Components of the second order Eshelby tensor for spheroidal inhomogeneities

We consider here the case of a spheroidal inhomogeneity. We put a 1 = a 2 = b and a 3 = a. By ǫ = a/b we define the aspect ratio of the spheroid. Expressions for I ij and T ijk are: (71)

I

a = T 1

 1 (a) + T 2 (a) + T 3 (a) + T 6 (a) + T 7 (a) + T 11 (a) + T 15 (a) +T 16 (a) + T 21 (a) + T 26 (a) + T 31 (a).

⊙ 3 T 7 T 8 T 9 T 15 T 7 T 7 T 8 7 T 8 T 9 T 10 T 10 T 11 10 T 11 T 12 T 13 T 13 T 14 13 T 14 T 15 Table 4 : 16 T 17 T 18 T 19 T T 20 T 21 T 22 20 T 21 T 22 T 23 T T 24 T 25 T 26 24 T 25 T 26 T 27 T T 28 T 29 T 30 28 T

 378915778789101011101112131314131415416171819202122202122232425262425262728293028 10 T 11 T 12 T 13 T 14 T The triple contraction between the T n for n = 7..15 ⊙ T 16 T 17 T 18 T 19 T 20 T 21 T 22 T 23 T 24 T 25 T 26 T 27 T 28 T 29 T 30 T 31 T T 16 T 17 T 18 T 19 0 29 T 30 T 31

T 3 (

 3 b) = a 3 T 3 (a) + a 4 T 4 (a), T 6 (b) = a 5 T 5 (a) + a 6 T 6 (a), T 7 (b) = a 7 T 7 (a) + a 8 T 8 (a) + a 9 T 9 (a), T 11 (b) = a 10 T 10 (a) + a 11 T 11 (a) + a 12 T 12 (a), T 15 (b) = a 13 T 13 (a) + a 14 T 14 (a) + a 15 T 15 (a), T 16 (b) = a 16 T 16 (a) + a 17 T 17 (a) + a 18 T 18 (a) + a 19 T 19 (a), T 21 (b) = a 20 T 20 (a) + a 21 T 21 (a) + a 22 T 22 (a) + a 23 T 23 (a), T 26 (b) = a 24 T 24 (a) + a 25 T 25 (a) + a 26 T 26 (a) + a 27 T 27 (a), T 31 (b) = a 28 T 28 (a) + a 29 T 29 (a) + a 30 T 30 (a) + a 31 T 31 (a) (36)

  B 3 }. We now introduce the tensor C defined by C = A⊙ 3 B; its components along (T 1 , ..., T 31 ) are denoted by c n for n = 1..31. The elements of the compact representation of tensor C are {c 1 , c 2 , C 1 , C 2 , C 3 } with:

2 , s 2 = I 13 4π b 2 , s 3 = 2 , s 5 = 2 , s 7 = 2 , s 9 =

 223252729 11 = I 22 = I 12 = π(3b 2 (η -1) + 2a 2 ) 2b 2 (a 2 -b 2 ) , I 13 = I 23 = 2π(1 -3η) a 2 -b 2 ,I 33 = 4π(3ηa 2 -b 2 ) 3a 2 (a 2 -b 2 ) , T 111 = T 112 = T 122 = T 222 = π(6a 2 b 2 (2 -3η) + 3b 4 (η -1) -4a 4 ) 12b 2 (a 2 -b 2 ) 2 , T 113 = T 123 = T 223 = π(3b 2 (η -1) + 12ηa 2 -2a 2 ) 2(a 2 -b 2 ) 2 , T 133 = T 233 = 2π(b 2 (5 -9η) -6a 2 η) 3(b 2 -a 2 ) 2 , T 333 = 4π(3a 2 (5η -1) -2b 2 )b 2 15a 2 (b 2 -a 2 )2 introduction of η allows to produce various equations which are valid for both the case of an oblate spheroid (a ≤ b) and a prolate spheroid (b ≥ a). Due to the symmetries of the problem, the second order Eshelby tensor is transversely isotropic and a representation of this tensor along the basis (T 1 , ..., T 31 ) is then possible: sn for n = 1..31 are defined by:s 1 = 3T 111 + 4(1 -ν 0 )I 11 4π(1 -ν 0 ) b T 113 + 2(1 -ν 0 )I 13 4π(1 -ν 0 ) a 2 , s 4 = T 113 + (1 -ν 0 )I 13 4π(1 -ν 0 ) b T 113 + (1 -ν 0 )I 13 2π(1 -ν 0 ) a 2 , s 6 = 2T 113 + (1 -ν 0 )(2I 11 + I 13 ) 4π(1 -ν 0 ) b T 113 + I 13 2π(1 -ν 0 ) a 2 , s 8 = 2T 113 + (1 -ν 0 )I 13 4π(1 -ν 0 ) b 3T 133 + 2ν 0 I 13 8π(1 -ν 0 ) a 2 , s 10 = 2T 113 + (1 + ν 0 )I 13 2π(1 -ν 0 ) a 2 , s 11 = 4T 113 + (1 -ν 0 )(4I 11 + I 13 ) 4π(1 -ν 0 ) b 2 , s 12 = 3T 133 + 2I 13 4π(1 -ν 0 ) a 2 , s 13 = 3[T 133 + 2ν 0 I 33 ] 4π(1 -ν 0 ) a 2 , s 14 = 3T 133 + 2I 13 -2ν 0 I 13 4π(1 -ν 0 ) b 2 , s 15 = 3[5T 333 + 2(2 -ν 0 )I 33 ] 8π(1 -ν 0 ) a 2 , s 16 = 3T 111 + (3 -2ν 0 )I 11 2π(1 -ν 0 ) b 2 , s 17 = 3T 111 + 2I 11 2π(1 -ν 0 ) b 2 , s 18 = T 113 + 2ν 0 I 11 4π(1 -ν 0 ) b 2 , s 19 = T 113 + (1 -ν 0 )I 13 2π(1 -ν 0 ) a 2 , s 20 = 3T 111 + I 11 4π(1 -ν 0 ) b 2 , s 21 = 3T 111 + 2I 11 4π(1 -ν 0 ) b 2 , s 22 = T 113 + 2ν 0 I 11 8π(1 -ν 0 ) b 2 , s 23 = T 113 4π(1 -ν 0 ) a 2 , s 24 = 2T 113 + ν 0 I 13 4π(1 -ν 0 ) b 2 , s 25 = T 113 + ν 0 I 13 2π(1 -ν 0 ) b 2 , s 26 = 3T 133 + 2(2 -ν 0 )I 13 8π(1 -ν 0 ) b 2 , s 27 = 3T 133 + 2(1 -ν 0 )I 13 4π(1 -ν 0 ) a 2 , s28= 4T 113 + (3 -ν 0 )I 13 8π(1 -ν 0 ) b 2 , s 29 = 2T 113 + (1 + ν 0 )I 13 4π(1 -ν 0 ) b 2 , s 30 = 3T 133 + 2I 13 8π(1 -ν 0 ) b 2 , s 31 = 3T 133 + (1 -ν 0 )(3I 33 + I 13 ) 4π(1 -ν 0 ) a 2 .

  1 , ..., T 31 ). The relations giving the coefficients a n for n = 1..31 as functions of the components A ijkpqr are given in appendix B. It must be emphasized that relations (60) can also be used to prove that tensors T n for n = 1..31 are linearly independent. Indeed, in (34), tensor A if the components A ijkpqr are null and, from relations (60), if a n = 0 for n = 1..31.The triple contraction between two elements of the basis T n are given in table 2 together with tables 3, 4 and 5. It can be observed that the set of tensors (T 1 , ..., T 31 ) is constituted of 5 groups which are independent for the composition ⊙ 3 . These groups are {T 1 }, {T 2 }, {T 3 , ..., T 6 }, {T 7 , ..., T 15 }, {T 16 , ..., T 31 }, and define sub-spaces of E 6 whose dimensions are respectively: 1, 1, 4, 9, 25. This irreducible basis is constituted of two kinds of elements:• the first kind of tensors are idempotent since they are invariants by the inner product "⊙ 3 ", i.e. T n ⊙ 3 T n = T n . These tensors are T 1 , T 2 , T 3 , T 6 , T 7 , T 11 , T 15 , T 16 , T 21 , T 26 , T 31 and correspond to the tensors that appear in the diagonal of table 1. • other tensors are nilpotent; the inner product "⊙ 3 " of these tensors by themselves gives the null tensor , ie. T n ⊙ 3 T n = 0. These tensors appear only out of the diagonal in table 1.

Table 2 :

 2 The triple contraction between the T n for n = 1..31

	⊙ 3 T 1 T 2 T 3	. . .	T 6 T 7	. . .	T 15 T 16	. . .	T 31
	T 1 T 1 0 0	. . .	0 0	. . .	0	0	. . .	0
	T 2 0 T 2 0	. . .	0 0	. . .	0	0	. . .	0
	T 3 0 0				0	. . .	0	0	. . .	0
	. . .	. . .	. . .		Table 3		. . .	. . .	. . .	. . .	. . .	. . .
	T 6 0 0				0	. . .	0	0	. . .	0
	T 7 0 0 0	. . .	0				0	. . .	0
	. . .	. . .	. . .	. . .	. . .	. . .		Table 4		. . .	. . .	. . .
	T 15 0 0 0	. . .	0				0	. . .	0
	T 16 0 0 0	. . .	0 0	. . .	0
	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .		Table 5
	T 31 0 0 0	. . .	0 0	. . .	0

Table 5 :

 5 The triple contraction between the T n for n = 16..31 By using the decomposition (35) for the third order tensor b = A ⊙ 3 a, we obtain:

  17 a 18 a 19 a 20 a 21 a 22 a 23 a 24 a 25 a 26 a 27 a 28 a 29 a 30 a 31

				a 16 a 
				
	   	, A 3 =	    	    	(38)
				

8 a 9 a 10 a 11 a 12 a 13 a 14 a 15

  5A ijkpqr π ij πpqπ kr /24 + A ijkpqr π ik πpqπ jr /12 -A ijkpqr π ik π jq πpr/3 + A ijkpqr π ir π kq πpr/3, a 2 = A ijkpqr π jq π kr n i np -A ijkpqr π jr π qk n i np,a 3 = A ijkpqr π ip π jq n k nr/2 -A ijkpqr π ij πpqn k nr/4, a 4 = A ijkpqr π iq π jr npn k /2 -A ijkpqr π ij πqrnpn k /4, a 5 = A ijkpqr π jp π kq n i nr -A ijkpqr π jk πpqn i nr/2, a 6 = A ijkpqr π jq π kr n i np/2 + A ijkpqr π jr π qk n i np/2 -A ijkpqr π jk πqrn i np/2, a 7 = A ijkpqr π ij πpqn k nr/2, a 8 = A ijkpqr π ij πqrnpn k /2, a 9 = A ijkpqr π ij npnqn k nr/2, a 10 = A ijkpqr π jk πpqn i nr, a 11 = A ijkpqr π jk πqrn i np, a 12 = A ijkpqr π ik npnqn j nr, a 13 = A ijkpqr πpqn i n j n k nr, a 14 = A ijkpqr πprn i n j n k nq, a 15 = A ijkpqr n i n j n k npnqnr, a 16 = A ijkpqr π ik π jq πpr/3 -A ijkpqr π ir π kq πpr/3, -A ijkpqr π ij πpqπ kr /6 + A ijkpqr π ip π jq π kr /6 + A ijkpqr π ik πpqπ jr )/6, a 17 = A ijkpqr π ik πpqπ jr /2 -A ijkpqr π ij πpqπ kr /4, a 18 = A ijkpqr π ik π jr npnq/2 -A ijkpqr π ij π kr npnq/4, a 19 = A ijkpqr π ik π jq npnr -A ijkpqr π ij π qk npnr/2, a 20 = A ijkpqr π ij πpqπ kr /4 -A ijkpqr π ip π jq π kr /4 -A ijkpqr π ik πpqπ jr /4 + A ijkpqr π ir π kq πpr/2, a 21 = 3A ijkpqr π ij πpqπ kr /8 -A ijkpqr π ik πpqπ jr )/4, a 22 = 3A ijkpqr π ij π kr npnq/8 -A ijkpqr π ik π jr npnq/4, a 23 = 3A ijkpqr π ij π qk npnr/4 -A ijkpqr π ik π jq npnr/2, a 24 = A ijkpqr π pk πqrn i n j /2, a 25 = A ijkpqr πpqπ kr n i n j /2, a 26 = A ijkpqr π kr n i n j npnq/2, a 27 = A ijkpqr π pk n i n j nqnr, a 28 = A ijkpqr π jp πrqn i n k /2, a 29 = A ijkpqr π jr πpqn i n k /2, a 30 = A ijkpqr π ir npnqn k n j /2, a 31 = A ijkpqr π ip n j nqn k nr.

	(60)

A Demonstration of the completeness and irreducibility of the T n (a)

We aim at proving that the third order tensors T n (a), defined by relations [START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF], are linearly independent. Consider then the following third order tensor c defined by:

Replacing in relations (36) tensor b by c and coefficients an by ϕn, the condition c = 0 leads to: (

Using the defintions for the T n (a) (see relations [START_REF] Walpole | Elastic behavior of composite materials : theoretical foundations[END_REF]), it can be easily show that the above relations also read:

Quantities α, β, γ, p, q, r, s, u, v s , v a , w being linearly independents, the above linear system is verified iff ϕn = 0 for n = 1..31.

C Derivation of the components of the sixth order Eshelby tensor

The expressions of the harmonic potential Φ i and of the biharmonic potential Ψ i can be found in Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] (see equations 12.13 and 12.14).They read:

For an interior point (within the inhomogeneity), Φr and Ψr are polynomial functions respectively of rank 3 and 5. The sixth order Eshelby tensor is then given by:

where the V i..j are given, for an interior point, by:

and the I i..j integrals are defined by:

The derivatives of V i..j read (see Mura 1987, eq. 11.40.3):

It can then be shown that:

where the T ijk are related to the I i..j integrals by:

The T ijk are invariant by any permutation of the indices i, j, k. An integral expression for the T ijk is also given in equation (48). It is obtained by replacing the I i..j in (67) by their expression given by (64). Finally, from (62) and (66), it is easy to obtain expression (47).