Wouter Kuijper

Victor Ermolaev

Olivier Devillers

Celestial Walk: A Terminating Oblivious Walk for Convex Subdivisions

Keywords:

We present a new oblivious walking strategy for convex subdivisions. Our walk is faster than the straight walk and more generally applicable than the visiblity walk. To prove termination of our walk we use a novel monotonically decreasing distance measure.

Point location in a convex subdivision is a classical problem of computational geometry for which several data structures have been designed with good complexities in the worst case [START_REF] Kirkpatrick | Optimal search in planar subdivisions[END_REF][START_REF] Chazelle | Fractional cascading: I. a data structuring technique[END_REF][START_REF] Franco | Planar point location revisited[END_REF]. These intricate solutions are often unused in favor of simpler algorithms based on traversal of the planar subdivisions using neighborhood relations between faces, also known as walking algorithms [START_REF] Bose | Online routing in convex subdivisions[END_REF][START_REF] Bose | Online routing in triangulations[END_REF][START_REF] Devillers | Walking in a Triangulation[END_REF]. These walking algorithms can also be used as a building block in randomized data structures for point location [START_REF] Ernst P Mücke | Fast randomized point location without preprocessing in two-and three-dimensional Delaunay triangulations[END_REF][START_REF] Devillers | The Delaunay Hierarchy[END_REF].

Amongst convex subdivisions, Delaunay triangulations received a lot of attention because of their practical importance. For Delaunay triangulations, essentially two walking strategies are used: the straight walk and the visibility walk [START_REF] Devillers | Walking in a Triangulation[END_REF]. The straight walk visits all faces crossed by a line segment between a known face and the query point, while the visibility walk goes from a face to another if the query point is on the side of the new face with respect to the supporting line of the edge common to the two faces (cf. Figure 1).

The straight walk trivially terminates in the face containing the query point and generalizes to any planar subdivision but with the inconvenience of not being oblivious: the starting face of the walk must be remembered during the whole walk.

The visibility walk is oblivious, but proving its termination requires the use of particular properties of the Delaunay triangulation, and actually the visibility walk may loop in other subdivisions [START_REF] Devillers | Walking in a Triangulation[END_REF] (cf. Figure 2).

Regarding performance, both walks may visit all triangles in the worst case and visit O( n) triangles when the points are evenly distributed [START_REF] Devroye | Expected time analysis for delaunay point location[END_REF][START_REF] Devillers | The worst visibility walk in a random Delaunay triangulation is O( n)[END_REF]. From a practical point of view, the visibility walk is simpler to implement and a bit faster in practice because it uses less predicates to walk from a triangle to its neighbor [START_REF] Devillers | Delaunay triangulations, theory vs practice[END_REF].

Contribution

We propose celestial distance 1 as a new way to measure the proximity between an edge of the subdivision and a query point. This distance measure allows to design new walking strategies 1 The name celestial distance refers to the practice of celestial navigation where angular distances between the celestial bodies and the horizon are used for navigation at sea.

RR n°9099

and prove their termination. By design these strategies are oblivious: only the current edge of the current face determines the edges by which the walk progresses to the next face.

Our main contribution is the celestial walk which is a refinement of the well-known visibility walk that has the additional advantage of terminating not just on Delaunay triangulations but on arbitrary convex subdivisions. This is particularly useful for constrained and/or incremental meshing where the conditions necessary for termination of the visibility walk can be locally and/or temporarily violated.

Another important feature of the celestial walk is that, like the visibility walk, it uses only orientation predicates to navigate the mesh. In practice, checking an orientation predicate reduces to computing the sign of a second degree polynomial, the degree of such polynomial being a relevant measure of the predicate complexity [START_REF] Liotta | Robust proximity queries: An illustration of degree-driven algorithm design[END_REF][START_REF] Boissonnat | Robust plane sweep for intersecting segments[END_REF]. As a consequence, it becomes relatively straightforward to implement the celestial walk in an efficient and robust manner.

Pre-requisites

Let G = (V, E, F ) be a planar straight line graph (PSLG) consisting of a set of vertices V , a set of half-edges E, and a set of faces F .

We abstract away the borders of G by assuming that it tiles the entire real plane. At the same time we rule out dense tessellations by assuming G is locally finite meaning the number of vertices (edges, faces) intersecting a given, bounded area is always finite.

We assume G is given in half-edge representation. In particular we assume the following atomic functions [START_REF] Mark De Berg | Computational Geometry: Algorithms and Applications[END_REF]:

origin ∶ E → V which maps every half-edge to its start vertex, target ∶ E → V which maps every half-edge to its end vertex, edge ∶ F → E which maps every face to some edge on its perimeter, face ∶ E → F which maps every half-edge to its corresponding (left-hand-side) face, next ∶ E → E which maps every half-edge to its successor half-edge in the counter-clockwise winding order of the face perimeter, twin ∶ E → E which maps every half-edge to its twin half-edge running in the opposite direction, i.e.: twin(twin(e)) = e, origin(e) = target(twin(e)) and vice versa.

Inria

The point location problem in PSLGs can now be formulated as follows: given some goal location p ∈ R 2 and an initial half-edge e init ∈ E, find some goal edge e goal ∈ E such that p ∈ face(e goal ) using only next(⋅) and twin(⋅) to get from one half-edge to the next, i.e.: there must exist a finite path e init = e 0 . . . e n = e goal such that for all 0 < i ≤ n it holds e i = next(e i-1 ) or e i = twin(e i-1 ).

Celestial Distance

One problem that we encounter when we try to use Euclidean distance as a measure of progress for a walking algorithm is the fact that Euclidean distance is not always strictly decreasing for every step in the walk. The latter means that it is not possible to prove termination using Euclidean distance alone. For this reason we define the following augmented distance measure on (half-)edges.

For a given point p ∈ R 2 and a half-edge e ∈ E we define the celestial distance of e to p as a pair cd(e, p) = [d, α] where d is the length of the line segment from p to the closest point on e and α is the wide angle (α ≥ π

2 ) between e and this line segment or 0 in case d = 0. We now define a lexicographic order on celestial distances as follows:

[d, α] < [d ′ , α ′ ] iff d < d ′ ∨ (d = d ′ ∧ α < α ′ )
We illustrate this new distance measure in Figure 3. The figure shows a convex face and the resulting partition of the plane obtained by grouping together points based on their closest edge. This leads to two superimposed Voronoi partitions: one based on classical Euclidean distance and one based on celestial distance. The partition based on celestial distance is a proper refinement of the partition based on Euclidean distance in the sense that points which were ambiguous under Euclidean distance are now partitionable under celestial distance. Note that the points in the corner-cones radiating outward from the vertices are all ambiguous under Euclidean distance (since their closest point on the polygon is the corner vertex which is shared by two edges) yet under celestial distance they can be further partitioned. In particular, for the partition based on celestial distance the ambiguous corner-cones are split according to the angular bisectors.

Abstract Walk

With all the pre-requisites and our celestial distance measure in place we are in a position to present the walking algorithm properly. We first present an abstract version of the algorithm and give a correctness proof for this abstract version.

In the next section, we will give an efficient, concrete instantiation of this abstract version where the computation of the celestial distances will be completely implicit. However, the correctness of the final version will rest on the correctness proof of the abstract version as given in this section.

The algorithm is rather simple: given a starting edge e and a query point p to the left of e, we select an edge of the face(e) that improves the distance to p. The abstract walk is formalized in Algorithm 1.

In lines 1-3 we bootstrap the invariant that the query point is always to the left of the current edge. In line 4 we bootstrap the current set of successor candidates. In line 5 we enter the main loop. The loop invariant ensures that the loop will terminate as soon as there are no more suitable successor edges that have lower celestial distance to the query point. In the proof below we will see how this condition is sufficient to ensure that, at termination, it holds p ∈ face(e). In closest point inside an edge of the polygon closest point is a vertex of the polygon angular bisector Proof. For termination note that, due to local finiteness of the mesh, and monotonicity of our distance measure, an infinitely descending chain of celestial distances is ruled out. It remains to prove that after termination it holds p ∈ face(e). For this it would suffice to show that for a half-edge e and goal location p such that p is to the left-of e but not in face(e) there always exists another e ′ in the perimeter of face(e) such that p is strictly on the right of e ′ and cd(e ′ , p) < cd(e, p). So let e p be the point on e closest to p and let r be the ray from p to e p . We Algorithm 1 Abstract Walk e and r (cf. Figure 4). Now it holds that: wideangle(e, r) = angle(e, e ′ ) + wideangle(e ′ , r) which implies wideangle(e ′ , r) < wideangle(e, r). This establishes cd(e ′ , p) < cd(e, p) and we conclude the algorithm would progress to twin(e ′ ).

Celestial Walk

In the previous section, the successor of an edge in the walk can be the twin of any edge of the current face with a smaller distance to the query point. In this section, we explain a way to actually select such an edge using only orientation checks.

The main problems that we are solving in this section are the facts that Algorithm 1 is nondeterministic and relies on the explicit computation of celestial distances. Both these properties make it less immediately applicable. In this section we therefore develop a derived algorithm that works for convex PSLGs. As we shall see, the additional assumption of convexity allows us a significantly more efficient walk.

So let us first consider the problem of determining a successor edge that has lower celestial distance than our current edge without having to explicitly compute these distances.

As an example of a convex face consider once more Figure 3. If we assume the query point is outside the face this leaves two possibilities. First, the query point may be located inside one of the orthogonal slabs (indicated in blue in Figure 3) which means the closest point coincides with the orthogonal projection of the query point on the edge. Second, the query point may be located in one of the intermediate corner-cones separating the orthogonal slabs which means the closest point coincides with the corner vertex.

The latter argument gives us a basic refinement of Algorithm 1 for the convex case: check if the query point is in one of the orthogonal slabs, if so, pick the corresponding edge, else the query point must be inside some corner-cone in-between two orthogonal slabs, in this case we can pick either edge (unless one of them is our current edge in which case we are forced to pick the other alternative in order to make progress).

To do even better we can resolve the remaining non-determinism by using the angle as a tie-breaker to determine which of the two candidate edges is best (i.e.: which edge is most favorably oriented towards the query point). Theoretically the best tie breaker is the angular bisector (indicated in purple in Figure 3). The only problem is that the explicit computation and representation of the exact angular bisector is at least as hard as the explicit computation and representation of celestial distances.

Fortunately we can avoid the explicit computation of angular bisectors as well. In particular we make the following case distinction. If the corner vertex is not obtuse (more precisely the face has an acute or right internal angle at that corner vertex) it holds that the extensions of both edges lie inside the corner-cone (i.e.: the leftmost corner in Figure 3). This means that we may use either edge itself as a crude approximation to the angular bisector (this is precisely what the visibility walk does in all cases). On the other hand, if the corner vertex is obtuse it holds that the normal to the base of the triangle spanned by the two candidate edges forms a suitable approximation to the angular bisector (because it lies properly inside the corner-cone). In Figure 5 we illustrate this construction. Note that, in the limit, as the internal angle approaches 180 degrees, or, alternatively, the ratio between the lengths of the two edges approaches 1, the base normal converges to the exact angular bisector.

The latter considerations lead us to define approx_bisector(e ′ , e ′′ ) for some pair of consecutive edges e

′ , e ′′ such that e ′′ = next(e ′ ) and angle(e ′ , e ′′ ) is obtuse. In particular, we let approx_bisector(e ′ , e ′′ ) denote the line from target(e ′ ) -which is equal by definition to origin(e ′′ ) -in the direction orthogonal to, and to the right of, the internal baseline segment that connects origin(e ′ ) and target(e ′′ ). For an illustration of this construction cf. Figure 5. We also define Inria obtuse(e ′ , e ′′ ) as a predicate that is true iff target(e ′′ ) is strictly right of the line from target(e ′ ) -which we assumed is equal to origin(e ′′ ) -in the direction orthogonal to, and to the left of, e ′ . For an illustration of this construction cf. Figure 6. With these additional definitions Algorithm 2 formalizes the celestial walk (see Figure 7). end if 16: end while 17: return face(e)

Algorithm 2 Celestial Walk

In lines 1-3 we bootstrap the invariant that the query point is always to the left of the current edge. In line 4 we initialize a perimeter edge variable used for iterating over the outline of the current face. In line 5 we enter the outer loop. As can be seen by inspecting the code-paths inside the loop body, the loop invariant ensures that the loop will terminate iff the query point is inside the current face. In line 6 we check if the query point lies strictly to the right of the perimeter edge. If this is the case we know, due to the face convexity, that the query point is outside the face. From that moment on the only goal is to find an edge that allows us to walk to a neighboring face. Therefore, in line 7 we introduce a second perimeter edge variable that will be used to iterate over pairs (e ′ , e ′′ ) of consecutive edges rather than singletons. In line 8 we state the negation of the condition that we are looking for: either the internal angle of the face is acute2 or, otherwise, the query point falls strictly to the right of the approximate bisector. In line 9 we shift the pair of edges along the face perimeter until said condition is reached. In line 11 we have exited the inner loop so we know that e ′ contains the edge that is "sufficiently closest" to the query point for us to continue the walk to the neighboring face. Hence we re-assign e and e

′ and drop back into the outer loop as though we would have restarted the algorithm on this new edge. The latter can be seen as a form of tail-recursion. Indeed, a recursive formulation is trivially possible by transforming the entire procedure into a function and substituting a recursive call at line 12. We have not opted for the recursive formulation because it would be slightly less efficient. The efficiency loss would be due to the bootstrapping check in lines 1-3 which, in the recursion, would be superfluous.

Conclusion

We have shown how our celestial distance measure allows to design novel walking strategies. In particular we improved the applicability of the popular visibility walk by refining it into what we call the celestial walk. Our walk is almost as simple to implement yet is guaranteed to terminate on arbitrary convex subdivisions.

If we assume that all obtuseness checks are pre-computed or memoized (at the expense of 1 or 2 bits per half-edge respectively) we can make the following observations w.r.t. the expected complexity of our walk.

The worst case for our walk materializes when all angles in the mesh are obtuse, i.e. in a regular hexagonal mesh, in this case we require 2 orientation tests per visited half-edge (the same number as the straight walk).

For triangulations a simple counting argument proves that there is at most a fraction of 1 3 of obtuse angles. This gives us a worst case rough estimate of 4 3 orientation tests per visited half-edge (better than the straight walk and only slightly worse than the visibility walk). In practice, we expect better performance since the fraction of obtuse angles will be much less than 1 3 in triangulations that exhibit some regularity.

For meshes that are guaranteed to be triangular, minor optimizations are possible by unrolling the loops of the general algorithm. Furthermore, it is a good heuristic to minimize the number of visited triangles by alternating clockwise and counterclockwise winding order in order to balance out the systematic errors in approximating the angular bisectors.

Implementation

A self contained prototype implementation and interactive demo is available [START_REF] Kuijper | Celestial walk demo[END_REF]. The implementation was done in type-script and heavily optimizes for readability and demonstratability rather than production level performance. The interactive demo part allows to compare the celestial walk to the visibility walk and the straight walk. In particular, it should allow the interested reader to get a feeling for the behavior and the relative strengths and weaknesses of these three algorithms.

Future Work

An obvious next question to consider is whether the walk can be generalized to 3 dimensions. We have some reason to believe this is the case.

In 2D as we progress according to the celestial walk we either move forward in the direction of the query point or orient ourselves towards it. Similarly, in 3 dimensions we expect to either progress in Euclidian sense or roll in the direction of the query point. This intuition leads to developing a suitable non-increasing distance measure quantifying this relationship. In particular to the walking problem in tetrahedralization we find that a combination of Euclidian distance as the first progression criterion and two successive angles are sufficient to guarantee the progress towards the tetrahedron containing the query point [START_REF] Kuijper | Zig-zagging in a triangulation[END_REF].

Finally, we would also be interested in obtaining experimental results regarding the trade-off that exists between steering (i.e.: choosing a next edge that is favorably oriented to the query point) and driving (i.e.: choosing a next edge quickly based on easily computed criteria).

Figure 1 :

 1 Figure 1: Straight (left) and visibility (right) walks

Figure 2 :

 2 Figure 2: Visibility walk may loop

Figure 3 :

 3 Figure 3: Voronoi diagrams of a polygon. In blue the segment Voronoi diagram of the open edges and the vertices, in purple the celestial Voronoi diagram of the edges.

Figure 5 :′′Figure 6 :

 56 Figure 5: Construction for approximating angular bisector.

Figure 7 :

 7 Figure 7: Celestial walk (obtuse angles marked by orange dots).

  Illustration for limit case of Theorem 1. make a case distinction on r. First assume that r intersects the interior of face(e). In this case it must hold that r intersects the boundary of face(e) at least one more time in another edge that is closer to p than the current edge for the Euclidean distance, and thus for our celestial distance (otherwise p would lie in the current face). Next assume that r does not intersect the interior of face(e), In this case it must hold that e p = origin(e) or e p = target(e). Assume, w.l.o.g., that e p = target(e). In this case we claim that it holds that e It then immediately follows that the closest point on e ′ to p must be the only point that e and e

				p
		e	′	r
		e
		Figure 4:
	2:	e ← twin(e)
	3: end if
	4: E 5: while E ′ ← {e} ′ ≠ ∅ do
	6: 7:	e ← select(E E ′ ← {twin(e ′ )
				Inria

1: if p strictly right of e then ′ ) | e ′ ∈ face(e) ∧ cd(e ′ , p) < cd(e, p) ∧ p strictly right of e ′ } 8: end while 9: return face(e) ′ = next(e) has strictly smaller celestial distance to p. Assume, w.l.o.g., that e ′ provides no improvement over e w.r.t. Euclidean distance to p. ′ share which, by assumption, is e p = target(e) = origin(e ′ ). Since r does not intersect the boundary of face(e) in any point other than e p it must hold that e ′ is in-between

RR n°9099

Note that the obtuseness check in line 8 can potentially be memorized by the mesh at the expense of one bit per half-edge.Inria

rue du Jardin Botanique CS20101 54603 Villers-lès-Nancy Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Salah Eddine Saidi
email: salah-eddine.saidi@ifpen.fr

Nicolas Pernet
email: nicolas.pernet@ifpen.fr

Yves Sorel
email: yves.sorel@inria.fr

AUTOMATIC PARALLELIZATION OF MULTI-RATE FMI-BASED CO-SIMULATION ON MULTI-CORE

Keywords: parallel, co-simulation, multi-core, acyclic orientation, scheduling

Co-simulation refers to simulating a complex system using several coupled numerical models. Engineers define the rate of data exchange between the models by setting communication steps. FMI is a standardized interface which easily allows coupling and co-simulation of numerical models. The RCOSIM approach allows the parallelization on multi-core processors of co-simulations using the FMI standard. In this paper, we tackle the limitations of this approach. First, we extend the co-simulation to multi-rate, i.e. with different communication steps. We present graph transformation rules and an algorithm that allow the execution of each model at its respective rate while ensuring correct data exchange between models. Second, we present an acyclic orientation heuristic for handling mutual exclusion constraints between operations that belong to the same model due to the non-thread-safe implementation of FMI. We evaluate the obtained speedup on a multi-core processor and the effect on the accuracy of the numerical results.

INTRODUCTION

The recent advancement in merging different technologies and engineering domains has led to the emergence of the field of Cyber-Physical Systems (CPS). In such systems, embedded computers interact with, and control physical processes. Because of the heterogeneity of the involved components (multi-physics, sensors, actuators, embedded computers), CPS may feature very complex architectures and designs, ever raising the need for time, cost, and effort-effective approaches for building robust and reliable systems. Numerical simulation has proven successful in responding to this need, and is therefore increasingly considered to be an indisputable step in design verification and validation. For complex CPS, experts of different engineering disciplines may be involved in the design process by developing models for different parts of the system. In a later stage, the developed models need to be integrated into one simulation environment and coupled together in order to perform simulation at the system level. This approach is called co-simulation. Each model is assigned an integration step, which in some cases is driven by the dynamics of the modeled system and the control objective, and exchanges data with the other models according to a communication step which can be equal or larger than its integration step. Enabling co-simulation of heterogeneous models requires using adequate tools and methods. In this scope, the Functional Mock-up Interface (FMI) [START_REF] Blochwitz | Functional Mockup Interface 2.0: The Standard for Tool Independent Exchange of Simulation Models[END_REF] was developed with the goal of facilitating the co-simulation and the exchange of numerical models. FMI defines a standardized interface that can be implemented by modeling tools in order to create models that can be connected with other FMI models. A model that is developed using FMI is exported as a Functional Mock-up Unit (FMU) which is a package that encapsulates an XML file containing, among other data, the definitions of the model's variables, and a library defining the equations of the model as C functions.

Co-simulation is an alternative approach to monolithic simulation where a complex system is modeled as a whole using differential equations and simulated by numerically solving these equations. In a CPS, the controlled physical process constitutes a multi-physics system and is modeled in the continuous-time domain using (hybrid) Ordinary Differential Equations (ODEs). Because they are aimed to be implemented on embedded computers, numerical laws that control the physical process are modeled in the discretetime domain. All of these features add to the complexity of the numerical models. Co-simulation has a number of advantages over the monolithic approach. It allows modeling each part of the system using the most appropriate tool instead of using a single modeling software. Also, it allows a better intervention of the experts of different fields at the subsystem level. Furthermore, co-simulation facilitates the upgrade, the reuse, and the exchange of models. In co-simulation, the equations of each FMU are integrated using a solver separately. FMUs exchange data by updating their inputs and outputs according to their communication steps. It is worth noting that in this paper, the term co-simulation is used to refer to the simulation of FMUs generated from FMI for Model Exchange, FMI for Co-Simulation, or a combination of both types. In our approach, when importing Model Exchange FMUs, a dedicated solver is assigned for each FMU. For connected FMUs with different communication steps, the master algorithm can provide extrapolation techniques to produce the missing data. In the following, we consider FMI co-simulation without step rejection.

Figure 1 shows the evolution of time and data exchange between three FMUs. The vertical double arrows represent data exchange between the FMUs. FMU B must be assigned a communication step H B equal to the smaller communication step among H A and H C since it exchanges data with both FMU A and FMU B. Co-simulation requires domain-specific information for each involved FMU. Such information, which is beyond the scope of FMI, can be provided by domain experts, for example by using an approach like the one propped by [START_REF] Sirin | A Model Identity Card to Support Simulation Model Development Process in a Collaborative Multidisciplinary Design Environment[END_REF]. The communication steps of FMUs can be shared as part of this information. Usually, assembling FMUs results in a heavy co-simulation, requiring high processing power. Increasing CPU performance through frequency scaling has reached some technological limits. As a consequence, during the last decade, parallelism has been by far the main way for increasing the performance of processors. In fact, the last years have witnessed a major shift among semiconductor manufacturers to building multi-core processors, i.e. integrating multiple processors into one chip allowing parallel processing on a single computer. Enabling parallel execution of heavy co-simulations on multi-core processors is keenly sought by the developers and the users of simulation tools. However, fulfilling this objective is not trivial and appropriate parallelization techniques need to be applied on co-simulation models in order to accelerate their execution on multi-core processors.

FMU A FMU B

In order to achieve (co-)simulation acceleration using multi-core execution, different approaches are possible and were already explored. From a user point of view, it is possible to modify the model design in order to prepare its multi-core execution, for example by using models that are developed using parallel computing libraries as in [START_REF] Gebremedhin | A Data-Parallel Algorithmic Modelica Extension for Efficient Execution on Multi-Core Platforms[END_REF]. Using parallel solvers as in [START_REF] Elmqvist | Parallel Model Execution on Many Cores[END_REF]) is another alternative. [START_REF] Galtier | FMI-based Distributed Multi-Simulation with DACCOSIM[END_REF] presented the DACCOSIM framework which allows multi-threaded distributed simulation of FMUs on multi-core processors as well as clusters of computers. However, the parallelization achieved through this approach is not automatic as the user has to define the distribution of the FMUs on the architecture. A more detailed discussion of related work is given in [START_REF] Ben Khaled | Fast Multi-core Co-simulation of Cyber-Physical Systems: Application to Internal Combustion Engines[END_REF].

In this paper, we address the problem from a co-simulation tool provider point of view. In such a tool, the user connects different FMUs, embedding solvers or using solvers provided by the co-simulation tool.

In this case, it is not possible to change the models, the solvers, nor the modeling tools. The Refined CO-SIMulation (RCOSIM) [START_REF] Ben Khaled | Fast Multi-core Co-simulation of Cyber-Physical Systems: Application to Internal Combustion Engines[END_REF]) approach allows the parallelization on multi-core processors of co-simulations using the FMI standard, by generating an efficient multi-core FMI master. This approach is based on a representation of the co-simulation using a Directed Acyclic Graph (DAG). Each vertex of this graph represents an operation which is an atomic unit of execution consisting of one or several FMU C function calls. Each edge represents a precedence relation between two operations.

In this paper, we deal with the limitations of this approach that we highlighted in a recent work [START_REF] Saidi | Acceleration of FMU Co-Simulation On Multi-core Architectures[END_REF]. Although RCOSIM resulted in interesting co-simulation speedups, it has two limitations that have to be considered in the multi-core scheduling problem in order to obtain better performances. First, RCOSIM supports only mono-rate co-simulations, i.e. the different FMUs have to be assigned the same communication step. Indeed, because of this limitation, RCOSIM is unable to handle an important number of industrial applications containing multiple communication steps as provided by the engineers. Second, the functions of an FMU may not be thread-safe, i.e. they cannot be executed in parallel as they may share some resource (e.g. variables) that might be corrupted if two operations try to use it at the same time. Consequently, if two or more operations belonging to the same FMU are executed on different cores, a mechanism that ensures these operations are executed in disjoint time intervals must be set up. Previously, these mutual exclusion constraints were tackled in RCOSIM by executing the operations of the same FMU on the same core which restricts the exploitation of the parallelism.

We propose in this paper a solution for each of the aforementioned limitations by making extensions to RCOSIM. In order to allow handling multi-rate co-simulations, we present a graph transformation algorithm that is applied to the initial constructed graph. Then, mutual exclusion constraints are represented by adding, in the graph, non oriented edges between operations belonging to the same FMU. We present an acyclic orientation heuristic that is used to assign directions to the added non oriented edges. The proposed solutions are implemented in the xMOD multi-model integration software. We evaluate the obtained speedup on a multi-core processor and assess the effect on the accuracy of the numerical results.

The rest of the paper is organized as follows. The next section gives an overview of the RCOSIM method and the proposed extensions. In section 3, graph transformation rules and an algorithm for handling multirate co-simulation are presented. In section 4, an acyclic orientation heuristic is presented. This heuristic is meant to handle the mutual exclusion constraints between the operations that belong to the same FMU. Section 5 describes the allocation and scheduling heuristic. In section 6, we present the obtained results.

We give a brief description of the industrial use case that is used to assess the proposed method and then present the obtained speedup and the numerical error evaluation. Finally, we conclude with some remarks and perspectives for future work in section 7.

OVERVIEW OF RCOSIM

The proposed method is an extension of the RCOSIM approach. We briefly present RCOSIM in this section before giving a description of the proposed extensions to tackle its limitations.

The entry point is a user-specified set of interconnected FMUs as depicted in Figure 2a. The execution of each FMU is seen as computing a set of input operations (one operation for each of the inputs of the FMU), a set of output operations (one operation for each of the outputs of the FMU), and one state operation for updating the state variables of the FMU. An operation is defined by a number of FMU C function calls. An input (resp. output ) operation is executed by calling fmiSet (resp. fmiGet) function and a state operation is executed by calling SetTime, GetDerivatives, SetContinuousStates, etc., functions in the case of FMI for Model Exchange or DoStep function in the case of FMI for Co-Simulation. Thanks to FMI, it is additionally possible to access information about the internal structure of a model encapsulated in an FMU.

In particular, as shown in Figure 2b, FMI allows the identification of Direct Feedthrough (e.g. Y B1 ) and Non Direct Feedthrough (e.g. Y A1 ) outputs of an FMU and other information depending on the version of the standard. RCOSIM consists in using FMU information on input/output dependencies to build a graph with an increased granularity, and then exploiting the potential parallelism [START_REF] Saidi | Acceleration of FMU Co-Simulation On Multi-core Architectures[END_REF]) of this transformed graph by using a heuristic to build an off-line multi-core schedule. The co-simulation can be described by a DAG G I (V I , A I ) called the operation graph where each vertex o i ∈ V I represents one operation and each oriented edge (o i , o j ) ∈ A I , that we call an arc, represents a precedence relation between two operations. The operation graph is built by exploring the relations between the FMUs and between the operations of the same FMU. A vertex is created for each operation and arcs are then added between vertices if a data dependency exists between the corresponding operations. If FMI 1.0 is used, which does not give information about the dependencies between the state variables computation and the input and output variables computations, we must add edges between all input operations and the state operation of the same FMU. Furthermore, edges connect all output operations and the state operation of the same FMU because the computation at the simulation step k of an output must be performed with the same value of the state as for all the outputs belonging to the same FMU. Running the co-simulation consists in executing the graph repeatedly. A new execution of the graph cannot be started unless the previous one was totally finished. The operation graph corresponding to the FMUs of Figure 2 is shown in Figure 3. In order to accelerate the co-simulation, the approach proceeds in two main phases: first, the operation graph is constructed by exploring the inter and intra-FMU dependencies and then, the operations are allocated to the available cores in such a way to minimize the makespan of the graph. The makespan corresponds to the execution time of the whole graph.

Y A1 Y A2 Y A3 U A1 U A2 U B1 U B2 U B3 Y B1 Y B2 FMU A FMU B (
Y A1 U B3 Ẋ A U A2 U B1 U A1 U B2 Ẋ B Y A2 Y B2 Y B1 Y A3
o i ) = R -S(o i ) -C(o i ) -E(o i ). R = max o i ∈V I (E(o i ))
denotes the critical path of the graph. These notations are used through the different phases of the proposed method in this paper.

Sections 3 and 4 describe new transformations that we propose to apply to the operation graph for handling multi-rate FMUs and mutual exclusion constraints between the operations beloning to the same FMU, respectively. Finally, the allocation of the different operations to the cores is achieved using a multi-core scheduling heuristic.

TRANSFORMATION ALGORITHM FOR MULTI-RATE GRAPHS

This section presents an algorithm that transforms the initial multi-rate operation graph G I (V I , A I ) into a mono-rate operation graph G M (V M , A M ). The aim of this transformation is to ensure that each operation is executed according to the communication step of its respective FMU and also to maintain a correct data exchange between the different FMUs, whether they have different or identical communication steps. Similar algorithms have been used in the real-time scheduling literature [START_REF] Ramamritham | Allocation and Scheduling of Precedence-Related Periodic Tasks[END_REF].

We define the Hyper-Step (HS) as the least common multiple (lcm) of the communication steps of all the operations: HS = lcm(H(o 1 ), H(o 2 ), . . . , H(o n )) where n = |V I | is the number of operations in the initial graph. The Hyper-Step is the smallest interval of time for describing a repeatable pattern of all the operations. The transformation algorithm consists first of all in repeating each operation o i , r i times where r i is called the repetition factor of o i and r i = HS H(o i ) . Each repetition of the operation o i is called an occurrence of o i and corresponds to the execution of o i at a certain simulation step. We use a superscript to denote the number of each occurrence, for instance o s i denotes the s th occurrence of o i . Then, arcs are added between operations following the rules presented hereafter. Consider two operations o i , o j ∈ V I connected by an arc (o i , o j ) ∈ A I , then adding an arc (o s i , o u j ) to A M , depends on the simulation steps (time) at which o s i and o u j are executed. In other words, if k s and k u are the simulation steps associated with o s i and o u j respectively, then the inequality k s ≤ k u is a necessary condition to add the arc (o s i , o u j ) to A M . In addition, o u j is connected with the latest occurrence of o i that satisfies this condition, i.e. with the occurrence o s i such that s = max(0, 1, . . . , r i -1) : k s ≤ k u . In the case where H(o i ) = H(o j ) (and therefore r i = r j ), occurrences o s i and o u j which correspond to the same number, i.e. s = u, are connected by an arc. On the other hand, if H(o i ) ̸ = H(o j ), we distinguish between two types of dependencies: we call the arc (o i , o j ) ∈ A I a slow to fast (resp. fast to slow)

dependency if H(o i ) > H(o j ) (resp. H(o i ) < H(o j )).
For a slow to fast dependency (o i , o j ) ∈ A I , one occurrence of o i is executed while several occurrences of o j are executed. In this case, arcs are added between each occurrence o s i : s ∈ {0, 1, . . . , r i -1}, and the occurrence o u j such that:

u = s × H(o i ) H(o j ) .
We recall that for a slow to fast dependency, the master algorithm can preform extrapolation of the inputs of the receiving FMU. For a fast to slow dependency (o i , o j ) ∈ A I , arcs are added between each occurrence o s i , and the occurrence o u j : u ∈ {0, 1, . . . , r j -1} such that:

s = u × H(o j ) H(o i ) .
Arcs are added also between the occurrences of the same operation, i.e. (o s i , o s ′ i ) where s ∈ {0, 1, . . . , r i -2} and s ′ = s + 1. Finally, for each FMU, arcs are added between the s th occurrence of the state operation, where s ∈ {0, 1, . . . , r i -2}, and the (s + 1) th occurrences of the input and output operations. The multirate graph transformation is detailed in Algorithm 1. The algorithm traverses all the graph by applying the aforementioned rules in order to transform the graph and finally stops when all the nodes and the edges have been visited.

M (V M , A M ) the new mono-rate graph; V M ← / 0; A M ← / 0; foreach operation o i ∈ V I do
Compute the repetition factor of o i :

r i ← HS H(o i ) ; Repeat the operation o i : V M ← V M ∪ {o s i }, s ∈ {0, . . . , r i -1}; end foreach arc (o i , o j ) ∈ A I do if H(o i ) > H(o j ) ((o i , o j ) is a slow to fast dependency) then Compute u = s × H(o i ) H(o j )
and add the arc

(o s i , o u j ) to the new graph G M (V M , A M ): A M ← A M ∪ {(o s i , o u j )}, s ∈ {0, . . . , r i -1}; end else if H(o i ) < H(o j ) ((o i , o j ) is a fast to slow dependency) then Compute s = u × H(o i )
H(o j ) and add the arc

(o s i , o u j ) to the new graph G M (V M , A M ): A M ← A M ∪ {(o s i , o u j )}, u ∈ {0, . . . , r j -1}; end else Add the arc (o s i , o s j ) to the new graph G M (V M , A M ): A M ← A M ∪ {(o s i , o s j )}, s ∈ {0, . . . , r i -1}; end end foreach operation o i ∈ V do
Add an arc between successive occurrences of o i :

A M ← A M ∪ {(o s i , o s+1 i )}, s ∈ {0, . . . , r i -2}; end foreach operation o i ∈ V such that T (o i ) = state do Add arcs (o s i , o s+1 j ) to the new graph G M (V M , A M ) between o i and every operation o j such that M(o j ) = M(o i ) and T (o j ) ∈ {input, out put}: A M ← A M ∪ {(o s i , o s+1 j )}, s ∈ {0, . . . , r i -2}; end

ACYCLIC ORIENTATION HEURISTIC FOR HANDLING MUTUAL EXCLUSION CONSTRAINTS

We describe in this section a heuristic for handling mutual exclusion constraints between operations belonging to the same FMU. These constraints are induced from the non thread-safe implementation of FMI. Without loss of generality, the superscript which denotes the number of the occurrence of an operation is not used in the rest of the paper for the sake of clarity. Each vertex of the graph G M (V M , A M ) represents an operation that we refer to using the notation o i .

In addition to the precedence relations represented by arcs, mutual exclusion relations are represented by (non oriented) edges in the mixed graph G X (V X , A X , D X ), where V X is the set of vertices (operations), A X the set of arcs, and D X the set of edges. This mixed graph is obtained from the graph G M (V M , A M ) by adding edges between operations belonging to the same FMU. Note that V X = V M . A precedence relation between the pair of operations o i and o j is denoted by an arc (o i , o j ) ∈ A which describes an order of execution between o i and o j , whereas a mutual exclusion relation between the pair of vertices o p and o q is denoted by an edge [o p , o q ] ∈ D X which means that o p and o q have to be executed in disjoint time intervals but in either order ([o p , o q ] = [o q , o p ]). Choosing one order or another does not impact the numerical results of the cosimulation since these operations do not have data dependencies. Still, we have to ensure mutual exclusion between them due to the non-thread-safe implementation of FMI. The goal of the heuristic proposed here is to assign directions to the edges of D X in order to define an order of execution between each pair of operations joined by an edge with the aim of minimizing the length of the critical path of the resulting graph, and subject to the constraint that no cycle is generated in the obtained graph as the term acyclic orientation suggests. Indeed, the existence of a cycle in the graph prevents the definition of an execution order for the operations belonging to this cycle, and leads to a deadlock during the execution. The output of the heuristic is a graph

G F (V F , A F , / 0) = G F (V F , A F ).
The problem of acyclic orientation of mixed graphs has been studied and shown to be NP-Hard in [START_REF] Andreev | Branch and Bound Method for Mixed Graph Coloring and Scheduling[END_REF][START_REF] Sotskov | Scheduling Problems and Mixed Graph Colorings[END_REF][START_REF] Al-Anzi | Using Mixed Graph Coloring to Minimize Total Completion Time in Job Shop Scheduling[END_REF]. In these studies, only the case of operations with unit (or equal) execution times were considered and the acyclic orientation problem was interpreted as a graph coloring problem. In this paper, we consider operations with different execution times and use an approach which takes into account that only operations belonging to the same FMU are joined by edges.

The orientation of the edges of D X means applying a function α :

{[o i , o j ] ∈ D X } → {(o i , o j ), (o j , o i )}. Without loss of generality, consider that α([o i , o j ]) = (o i , o j ), then o i is added to Γ(o j ) and o j is added to Γ(o i ).
Consequently, the inequality E(o i ) ≤ S(o j ) holds and therefore o i and o j are executed in disjoint time intervals. The first step of the heuristic consists in building the mixed graph G X (V X , A X , D X ) by adding edges between operations belonging to the same FMU. The output of the heuristic is a topological ordering of the operations of each FMU. This can be interpreted for each FMU M as an assignment of integers φ : V F → {1, 2, . . . , m} defining the order of execution for these operations, where m is the number of operations of the FMU M . We call the integer assigned to an operation its rank.

Algorithm 2 details the acyclic orientation heuristic. After adding the edges to the graph, the heuristic iteratively assigns ranks to the operations. It keeps a list of already ordered operations l for each FMU M and at each main iteration it selects the operation which has the earliest start time S(o i ) to be added to the list of its FMU. Ties are broken by selecting the operation with the least flexibility F(o i ). We call the operation to be ordered at a given iteration, the pending operation. The the pending operation o i is assigned the ranks τ ∈ {x, x + 1, . . . , length(l) + 1} in ascending order, where x = max o j ∈Γ(o i ) (φ (o j ) + 1) (1 i f Γ(o i ) = / 0). For each rank τ, the assignment φ (o i ) = τ is made and the ranks assigned to all the already ordered operations

o i ′ ∈ l such that φ (o i ′ ) ≥ φ (o i ) are increased φ (o i ′ ) = φ (o i ′ ) + 1. Then, every edge [o i , o ′ i ] ∈ D X : o ′ i ∈
l is assigned a direction from the operation with the lower rank to the operation with the higher rank. At this point, the increase in R, the critical path of the graph, is evaluated. Next, the operations o i ′ ∈ l such that φ (o i ′ ) > φ (o i ) are reassigned their previous ranks φ (o i ′ ) = φ (o i ′ ) -1, and the pending operation is assigned the next rank φ (o i ) = φ (o i ) + 1. The increase in the critical path is evaluated again similarly. After repeating this for all the ranks τ ∈ {x, x + 1, . . . , length(l) + 1}, the pending operation is finally ordered at the rank which leads to the least increase in the critical path and edges

[o i , o ′ i ] ∈ D X : o ′ i ∈ l
are assigned directions as described before. The heuristic begins another iteration by selecting a new operation to be ordered. The heuristic assigns a rank to one operation at each iteration. Every operation is assigned a rank higher than the ranks of all its predecessors which guarantees that no cycle is generated. The heuristic finally stops when all the operations have been assigned ranks.

ALLOCATION AND SCHEDULING HEURISTIC

In order to achieve fast execution of the co-simulation on a multi-core processor, an efficient allocation and scheduling of the operation graph has to be performed. A multi-core scheduling heuristic is applied on the final graph G F (V F , A F ) obtained after the acyclic orientation phase. Because of the number of fine-grained operations, and since the execution times and the dependencies between the operations are known before predecessors that were allocated to different cores. The heuristic repeats this procedure and finally stops when all the operations have been scheduled.

EVALUATION

The proposed parallelization approach has been applied on an industrial use case. In this section, we give a description of the use case and then present the obtained results.

Tests have been performed on a computer running the Windows 7 operating system with 16 GB RAM and an Intel core i7 processor running 8 cores at 2.7 GH. Experiments have been carried out on a Spark Ignition (SI) RENAULT F4RT engine co-simulation. It is a four-cylinder in line Port Fuel Injector (PFI) engine in which the engine displacement is 2000 cm 3 . The air path is composed of a turbocharger with a mono-scroll turbine controlled by a waste-gate, an intake throttle and a downstream-compressor heat exchanger. This co-simulation is composed of 5 FMUs: one FMU for each cylinder and one FMU for the airpath. The FMUs were imported into xMOD using the FMI export features of the Dymola tool. This use-case leads to an initial graph containing over 100 operations. We refer to our proposed method as MUO-RCOSIM (for Multi-Rate Oriented RCOSIM). We compared the obtained results with two approaches: The first one is RCOSIM which is mono-rate and thus we had to use the same communication step for all the FMUs. We used a communication step of 20µs. The second one consists in using RCOSIM with the multi-rate graph transformation algorithm. We refer to it as MU-RCOSIM (for Multi-Rate RCOSIM). For MUO-RCOSIM and MU-RCOSIM we used the recommended configuration of the communication steps for this use case. For each cylinder, we used a communication step of 20µs. The communication step used for the airpath is 100µs. In fact, the airpath has slower dynamics than the cylinders and this configuration of the communication steps corresponds to the specification given by engine engineers. For each FMU, we used a Runge-Kutta 4 solver with a fixed integration step equal to the communication step. The graph of this use case is transformed by Algorithm 1 into a graph containing over 280 operations that are scheduled by the multi-core scheduling heuristic. The validation of the numerical results of the co-simulation using the proposed method is achieved through the comparison of the co-simulation outputs with reference outputs. Since it is not possible to solve the equations of the FMU analytically, the reference outputs are obtained by using RCOSIM which has been shown in [START_REF] Ben Khaled | Fast Multi-core Co-simulation of Cyber-Physical Systems: Application to Internal Combustion Engines[END_REF] to give a very good accuracy of the numerical results. Figure 5a shows the obtained results for the torque (an output of the airpath). We note that the results match with the reference, and the generated error is very small remaining within an acceptable bound (< 1%). Similar accuracy results were obtained for the different outputs of the co-simulation.

The speedup obtained using MUO-RCOSIM is compared with the speedups obtained using RCOSIM and MU-RCOSIM. The speedup was evaluated by running the co-simulation in xMOD. Execution times measurements were performed by getting the system time stamp at the beginning and at the end of the cosimulation. For a given run of the co-simulation, the speedup is computed by dividing the mono-core co-simulation execution time of RCOSIM by the co-simulation execution time of this run on a fixed number of cores. Figure 5b sums up the results. The same speedup is obtained using MUO-RCOSIM and MU-RCOSIM even when only 1 core is used. This speedup is obtained thanks to using the multi-rate configuration. More specifically, increasing the communication step of the airpath from 20µs to 100µs results in fewer calls to the solver leading to an acceleration in the execution of the co-simulation. By using multiple cores, speedups are obtained using both MUO-RCOSIM and MU-RCOSIM. Additionally, MUO-RCOSIM outperforms MU-RCOSIM with an improvement in the speedup of, approximately 30% when 2 cores are used, and approximately 10% when 4 cores are used. This improvement is obtained thanks to the acyclic orientation heuristic which defines an efficient order of execution for the operations of each FMU that are mutually exclusive. This defined order tends to allow the multi-core scheduling heuristic to better adapt the potential parallelism of the operation graph to the effective parallelism of the multi-core processor (number of cores) resulting in an improvement in the performance. MU-RCOSIM, on the other hand, uses the solution of RCOSIM which consists in simply allocating mutual exclusive operations to the same core introducing restrictions on the possible solutions of the multi-core scheduling heuristic. When using 8 cores, no further improvement is possible since the potential parallelism is fully exploited. Worse still, the overhead of the synchronization between the cores becomes counter-productive, which explains why the speedup with 8 cores is less than the speedup with 4 cores for all the approaches. The best performance is obtained using 5 cores with slight improvement compared to using 4 cores.

CONCLUSION

This paper addressed the problem of automatic parallelization of multi-rate co-simulations using the FMI standard. We proposed an extension of the RCOSIM approach. First, we presented a set of rules and an algorithm in order to allow the execution of multi-rate FMUs where each FMU is executed according to its communication step, and correct data exchange between FMUs is ensured. Second, we presented an acyclic orientation heuristic that is applied on the operation graph in order to define an execution order between the operations of the same FMU which are mutually exclusive due to the non-thread-safe implementation of FMI. This heuristic results in a graph which can be efficiently exploited by the multi-core scheduling heuristic leading to an improvement in the co-simulation performance. The obtained speedup and numerical accuracy results show the efficiency of our proposed method. In future work, we aim at comparing our solution with an exact scheduling algorithm and also an online scheduling approach. Also, we will test our proposed method on other industrial applications. In addition, we envision to extend MUO-RCOSIM to real-time multi-core scheduling in order to perform Hardware-in-the-Loop simulation whose main challenge is to map the real-time constraints on the different operations of the graph.

Figure 1 :

 1 Figure 1: Evolution of time and data exchange between three FMUs during co-simulation

  Figure 2: An example of inter and intra-FMU dependencies of two FMUs connected by the user

Figure 3 :

 3 Figure 3: Operation graph obtained from the FMUs of Figure 2

Figure 4 Figure 4 :

 44 Figure4shows the graph obtained by applying the multi-rate transformation algorithm on the graph of Figure3. In this example H B = 2 × H A , where H A and H B are the communication steps of FMUs A and B respectively.

Figure 5 :

 5 Figure 5: Speedup and numerical results of e-RCOSIM

Algorithm 2: Acyclic orientation heuristic Input: Operation graph G M (V M , A M ); Output: New operation graph G F (V F , A F ); Set G M (V M , A M ) the operation graph, G X (V X , A X , D X ) the mixed graph, and G F (V F , A F ) the final graph; V X ← V M ; A X ← A M ; D X ← / 0; V F ← V M ; A F ← A M ; Set Ω the set of all the operations:

. Break ties by selecting the operation with the least flexibility; Set M ← M(o i ) and l the list of the already ordered operations of FMU M ; Set σ ← ∞; (Initialize the increase in the critical path);

is more convenient to use an offline scheduling heuristic which has the advantage of introducing lower overhead than online scheduling heuristics. We use an offline scheduling heuristic similar to the one proposed in [START_REF] Grandpierre | Optimized Rapid Prototyping For Real-Time Embedded Heterogeneous multiprocessors[END_REF]) which is a fast greedy algorithm whose cost function corresponds well to our minimization objective. The heuristic considers the different timing parameters of each operation o i ∈ V F in order to compute a schedule that minimizes the makespan of the graph. The heuristic schedules the operations of the graph G F (V F , A F ) on the different cores iteratively and aims at minimizing a cost function which expresses the schedule pressure of an operation on a specific core. The schedule pressure is the difference between the makespan increase, by allocating this operation to this core, and the operation's flexibility. The heuristic updates the set of candidate operations to be scheduled at each iteration. An operation is added to the set of candidate operations if it has no predecessor or if all its predecessors have already been scheduled. For each candidate operation, the schedule pressure is computed on each core and the operation is allocated to its best core, the one that minimizes the pressure. Then, a list of candidate operation-best core pairs is obtained. Finally, the operation with the largest pressure on its best core is selected and scheduled. Synchronization operations are added between the scheduled operation and all its AUTHOR BIOGRAPHIES SALAH EDDINE SAIDI is a PhD student at Pierre et Marie Curie University and holds a Masters in Computer Engineering from the same university. As part of the PhD program he is a research scholar at IFP Energies nouvelles and INRIA. His research interests lie in co-simulation parallelization, multi-core, embedded, and real-time computing. His email address is salah-eddine.saidi@ifpen.fr.

NICOLAS PERNET is a Software Project Manager at IFP Energies nouvelles. He holds a PhD in Computer Science from Pierre et Marie Curie University. His research interests include real-time systems, multi-core programming paradigms, and co-simulation approaches. His email address is nicolas.pernet@ifpen.fr.

YVES SOREL is Research Director at Inria, the French National Institute for computer science and applied mathematics, and scientific leader of the team AOSTE-Paris. His main research topic concerns the analysis and the optimization of distributed real-time embedded systems. He is also the founder of the AAA methodology (Algorithm Architecture Adequation) and of SynDEx, a system level CAD software for distributed real-time embedded systems. His email address is yves.sorel@inria.fr.


Brice Le Borgne

Olivier De Sagazan

Samuel Crand

Emmanuel Jacques

Maxime Harnois

Conformal Electronics Wrapped Around Daily-life Objects Using Original Method: Water Transfer Printing

Keywords: Water Transfer Printing, Conformal Electronics, Curvilinear Electronics, Flexible electronics, Additive Transfer

The water transfer printing method is used to transfer patterned films on random threedimensional objects. This industrially viable technology has been demonstrated to intimately wrap metallic and polymeric films around different materials. This method avoids the use of rigid substrate during the transfer step. Patterns can be transferred to objects without folds even when holed, addressing a challenging issue in the field of conformal electronics. This technique enables to reach high films bending properties. This promising method, enables to integrate large-area films onto daily-life objects. Bent capacitive touchpad is fabricated showing the potential applications of this technology.

Most of electronic circuits are fabricated using rigid and brittle materials patterned on twodimensional (2D) planar substrates. Indeed, materials and fabrication methods, commonly used in microelectronics, ranging from lithography to materials etching or growth are usually incompatible with daily-life three-dimensional (3D) objects. However, countless new applications require curvilinear objects with built-in electronics systems.

Consequently, new materials and methods have to be developed to answer a crucial issue: how to integrate electronics which intimately wrap 3D surfaces which can be: soft (organics), curvilinear (daily-life objects) or stretchable (textiles). The past few years, answering this issue has gained a growing interest and many proofs of concept have been demonstrated including epidermal electronics, 1, 2, 3 medical devices, [START_REF] Kim | Dissolvable Films of Silk Fibroin for Ultrathin Conformal Bio-Integrated Electronics[END_REF][START_REF] Kim | Stretchable Silicon Nanoribbon Electronics for Skin Prosthesis[END_REF] curvilinear electro-optics, [START_REF] Song | Digital Cameras with Designs Inspired by the Arthropod Eye[END_REF][START_REF] Li | Mechanics of Curvilinear Electronics and Optoelectronics[END_REF] conformal displays, [START_REF] Sekitani | Stretchable Active-Matrix Organic Light-Emitting Diode Display Using Printable Elastic Conductors[END_REF][START_REF] Vosgueritchian | Stretchable LEDs: Light-Emitting Electronic Skin[END_REF] conformal photovoltaic devices, [START_REF] Kaltenbrunner | Ultrathin and Lightweight Organic Solar Cells with High Flexibility[END_REF] or wearable electronics. [START_REF] Rogel | Spontaneous Buckling of Multiaxially Flexible and Stretchable Interconnects Using PDMS/Fibrous Composite Substrates[END_REF][START_REF] Stoppa | Wearable Electronics and Smart Textiles: A Critical Review[END_REF][START_REF] Park | Graphene-Based Conformal Devices[END_REF] To date, most of conformal electronic devices are fabricated on 2D planar substrates which can be flexible or stretchable prior to the transfer step to the 3D object. In the case of flexible substrate, the thinner the substrate is, the more the devices conform to the object. [START_REF] Salvatore | Wafer-Scale Design of Lightweight and Transparent Electronics that Wraps Around Hairs[END_REF][START_REF] Feiner | Engineered Hybrid Cardiac Patches with Multifunctional Electronics for Online Monitoring and Regulation of Tissue Function[END_REF][START_REF] Fukuda | Fully-Printed High-Performance Organic Thin-Film Transistors and Circuitry on One-Micron-Thick Polymer Films[END_REF] However, wrapping capabilities improvement is at the expense of processing and handling convenience, especially for large-area applications. In the case of a stretchable substrate, its elongation capability allows intimate contact of devices with the object surface. However, electronic devices are often made of brittle materials which show moderate elongation capability. Consequently, elongation mismatches exist between the substrate and the materials which dramatically increase devices mechanical failure. Significant works have been performed to overcome this challenging issue including: new fabrication methods, [START_REF] Kim | Dissolvable Films of Silk Fibroin for Ultrathin Conformal Bio-Integrated Electronics[END_REF][START_REF] Vosgueritchian | Stretchable LEDs: Light-Emitting Electronic Skin[END_REF][START_REF] Rogel | Spontaneous Buckling of Multiaxially Flexible and Stretchable Interconnects Using PDMS/Fibrous Composite Substrates[END_REF][START_REF] Ko | Curvilinear Electronics Formed Using Silicon Membrane Circuits and Elastomeric Transfer Elements[END_REF] new mesh design, [START_REF] Zhang | Buckling in Serpentine Microstructures and Applications in Elastomer-Supported Ultra-Stretchable Electronics with High Areal Coverage[END_REF] or highly stretchable materials development. [START_REF] Sekitani | Stretchable Active-Matrix Organic Light-Emitting Diode Display Using Printable Elastic Conductors[END_REF] Importantly, using previously described techniques, conformal electronics wrapping have been demonstrated on curvilinear objects (concave or convex), even if folds occur systematically. However, wrapping holed objects remains challenging. Indeed, electronics cannot easily conformally wrap a ring because of the 2D substrate (see in the figure S1 of the supplementary file), even if it is water soluble. Indeed, in this particular case, folds will occur during the positioning step and remain after the transfer. Consequently, new transfer methods must be proposed in order to address the challenging issues previously mentioned.

In this work, the concept of water transfer printing is detailed. Besides, this method offers interesting bendability and wrapping capabilities. They are demonstrated for several kinds of 3D objects with various materials and morphological parameters such as roughness and shapes (from curvilinear to cubic). Moreover, interconnects printed on large-area substrate have been successfully transferred to hand-sized plastic object. Finally, metal/insulator/metal (MIM) capacitors and capacitive sensors that could potentially be used to detect finger position for touchpad applications have been transferred to curved surface. Note that, previous studies reported techniques close to the one presented here. Nevertheless, these works proposed to dissolve the substrate after the transfer step which does not prevent folds. [START_REF] Yeo | Multifunctional Epidermal Electronics Printed Directly Onto the Skin[END_REF][START_REF] Lee | Stretchable Graphene Transistors with Printed Dielectrics and Gate Electrodes[END_REF] An other uses copper substrate that must be etched by metal etchant that limits materials and methods commonly used in microelectronics. [START_REF] Zhang | Graphene Water Transfer Printing for 3D Surface[END_REF] Here, the process uses water-soluble films (PVA: polyvinyl(alcohol)) containing the active layers that are gently deposited on top of the water. In a few seconds, the PVA is dissolved and patterns float on water due to surface tension. Then, the 3D object is dipped into the water through the patterns. Then, layers wrap the surface because water pressure forces materials to spread on the object surface. However, a challenging issue, overcome in this work, relies on thin film patterning on PVA substrate. Indeed, processing on water-soluble substrate is not trivial. For instance, photolithography steps need chemical reagents that contain water in most of the cases. This study answers these technical issues and give several processing routes which are detailed in the experimental section (see supplementary file) and shown in figures S2 and S3. Thus, wide range of applications requirements, materials and methods, can be addressed to allow the fabrication of electronics directly on water-soluble films ranging from lithographical accuracy (i.e., small-area processing at high accuracy) to printing accuracy (i.e., large-area processing at medium accuracy). Indeed, the process depicted in figure S2 demonstrates the possibility to pattern evaporated thin film at micrometric scale on PVA substrate using conventional lift-off technique. Moreover, another process depicted in figure S3 shows that additive technologies such as inkjet-printing are a convenient way to achieve large-area processing on water-soluble films. Note that, the broad diversity of processes (figures S2 and S3) constitute a novel finding compared to the state of the art.

Figure 1 and Video S1 of supplementary file summarize the technological steps of water transfer printing technology. Figure 1a shows a 3D scheme of multilayers films (generic structure) patterned on PVA. Note that, the structures can be made of metal and/or polymer using photolithography and/or printing techniques over small or large-area. In order to illustrate the concept, the multilayered structure is composed of two layers only. The first layer is aluminum ribbons (150µm width and 1cm length; 150nm thick) ended by a square (2mm sidelength) designed in an 8 pointed regular star configuration. The second layer is a polymer designed in squared mesh (1mm side-length) using inkjet-printed Su8 (2000.5 MicroChem). In this work, Su8 has been chosen for its transparency but also because it can be easily printed, [START_REF] Robin | Epoxy Based Ink as Versatile Material for Inkjet-Printed Devices[END_REF] or patterned using conventional lithography on PVA substrate (water-free developer). Such interesting patterned structure will be discussed in the following section. Figure 1b and1e show the second process step. The PVA substrate is deposited on top of water poured into a dip tank.

Optical picture in figure 1f shows that the PVA substrate floats on the water. The dissolution of the PVA film occurs directly when the substrate is put into contact with water. As shown on the optical picture of figure 1g, the dissolution is firstly observed on the substrate edges and after less than 1 minute, the patterns float on top of the liquid interface (figure 1h). This mixture, composed of PVA dissolved in water, shows highly adhesive and viscous characteristics and remains on top of the water. Schematics in figure 1c and optical picture in figure 1i illustrate the dipping step. The 3D object is dipped through the floating patterns that intimately wrap its surface due to water pressure. The object is gently shacked into the water bath before its withdrawing (see figure 1d and1j). The object is then dried and residues of PVA that form a capping layer, can be removed by immersing a second time the object into the water. Then, the desired patterns firmly stick the object due to Van Der Waals forces. Note that, alignment accuracy between patterns and 3D object can be improved as shown in supplementary video S2 using 3-axis machine coupled with a visualization system. In this case, the user can proceed to alignment during the dissolution step at millimeter accuracy.

In the case described in figure 1, a structural polymer thin film is printed above the aluminum patterns. As shown in Figure S4, the polymeric mesh avoids the movement of aluminum during the PVA dissolution. Indeed, if mesh is not used, the aluminum patterns are not linked anymore after substrate dissolution. Radial forces then occur which drastically deform the final design (see figure S4b compared to figure S4a). Note that, the structural layer (polymer mesh) can also be functional leading to free-standing electronics. [START_REF] Feiner | Engineered Hybrid Cardiac Patches with Multifunctional Electronics for Online Monitoring and Regulation of Tissue Function[END_REF] In this case patterns can be fabricated at lithographic scale. Indeed, electrodes spaced from 50µm have been successfully transferred to curved glass, as shown in figure S5. However, this configuration does not allow to suppress folds on complex 3D objects such as angular shapes. Moreover, conformal wrapping remains challenging for holed objects.

To solve this challenging issue, water transfer printing technique can be improved. As shown in Figure S4c, rigid film guides can be placed on the edges of the PVA substrate in order to counter radial forces leading to a drastic decrease of patterns deformation. Using such improvement, aluminum patterns can remain unlinked and no additional layer acting as flexible substrate is necessary. Consequently, patterns can be transferred without folds, inside holes of complex 3D objects. Such assertions will be demonstrated in the next section and can be considered as another innovative part of this work.

As for other transfer printing techniques, hydrodipping enables the adhesion between the object and the patterns thanks to non-covalent force. In this case, surface energy is a critical experimental parameter which governs the strength of adhesion. [START_REF] Hur | Nanotransfer Printing by Use of Noncovalent Surface Forces: Applications to Thin-Film Transistors that Use Single-Walled Carbon Nanotube Networks and Semiconducting Polymers[END_REF] Consequently, in order to demonstrate the robustness of the water transfer printing method, aluminum ribbons have been transferred to 3D objects showing a broad range of surface energy. Standard glass (Pyrex), PMMA (poly(methyl methacrylate)) and Teflon have been chosen for their well-known surface energy values equal 68 mN/m, 23 38 mN/m 24 and 22 mN/m, [START_REF] Fu | Surface Composition and Surface Energy of Teflon Treated by Metal Plasma Immersion Ion Implantation[END_REF] respectively. Optical pictures in figures 2a, 2b and 2c show that aluminum ribbons are successfully transferred whatever the material wetting behavior of the 3D object is (water contact angle values are equal to 34°, 69° and 104°, for glass, PMMA and Teflon, respectively). Interestingly, aluminum ribbons can be successfully transferred to the curved shapes (concave and convex) of a glass beaker, the edge of a PMMA cube or the face of a Teflon cube. Moreover, the transfer process is reproducible and robust for these materials.

Figure 2d and 2e, show Scanning electron microscope (SEM) images of aluminum ribbon transferred on a 50 cents euro coin. Figure 2d shows that the aluminum ribbon is firmly anchored along the edge of a star. The maximum thickness of the star equals 70µm (see 3D profile figure S6a) and wrapping capability is demonstrated. However, inevitably, micro-folds occur when film coats an arbitrary 3D complex surface. A zoom on aluminum ribbon edge (figure 2e) shows that different behaviors can be distinguished when highly different morphologies are involved. The area marked by the number 1 shows a low and relatively homogeneous roughness (RMS equals to 0.8µm; see 3D profile in figure S6b). In that case, aluminum ribbon wraps intimately the object morphology. The area marked by the number 2 shows a steep step on the surface and the aluminum ribbon is locally released from the object.

Interestingly, the area marked by the number 3 shows a sharp ridge that induced extreme bending radius as low as several microns (see in 2D profiles in figure S6c,d,e,f andg). Note that, such interesting bending capabilities are allowed because no rigid substrate is used during the transfer step. Indeed, according to the strain theory, [START_REF] Gleskova | Failure Resistance of Amorphous Silicon Transistors under Extreme In-Plane Strain[END_REF] the strain level induced in the aluminum ribbon during the transfer is lower using water transfer printing than using conventional flexible substrate. Moreover, as shown in figure S7, interconnects have been transferred on human hair without mechanical and electrical failure.

Figure 2f shows the optical picture of aluminum ribbons that conformally wrapped the Teflon ring. The outside and inside diameters of Teflon ring equal 5cm and 1.4cm, respectively. The ring is a convenient shape to demonstration the advantages of water transfer printing technology. As shown in the figure S1 of the supplementary file, aluminum ribbons patterned on a flexible substrate cannot intimately wrap a holed surface such as a ring. Indeed, the conformal behavior suffers from the flexible substrate that covers the hole. Such drawback does not exists using the water transfer printing technique because the rigid substrate is dissolved during the transfer step. Note that, in this case, no polymeric mesh layer has been added to the star shaped pattern in order to avoid additional layer acting as flexible substrate. Figure 2g shows inkjet-printed silver (Ag) interconnects on A4 sized PVA substrate (21×29.7 cm). Water transfer printing has been successfully performed to wrap metallic wires around hand-sized plastic model car. Figures 2h and2i show optical pictures of the model car before and after the transfer step proving large-scale processing capability.

In order to demonstrate the water transfer printing technique potential applications, usual passive devices (MIM capacitors) have been successfully transferred to a bent shape (radius of curvature equals to 5 cm, Figure 3a) and to the edge (approximately 90°) of a plastic object (Figure 3b). These capacitors are composed of two 150nm-thick aluminum layers (top and bottom electrodes) and a 1.3 µm polymeric layer (dielectric) leading to a capacitance of about 34pF. Electrical characteristics after the transfer step are reported in Figure S8. Finally, a capacitive touchpad has been fabricated on a bent plastic object. The 3D object bending radius equals 5cm. Figure 3c shows the capacitive touchpad conformally wrapped around the object.

The sensing area (4cm×1cm) is defined by eight aligned electrodes (5mm×1cm) in order to detect finger sliding. V-shaped electrodes design is used to optimize the overlap between 1mm spaced electrodes. In order to respect the inter-electrodes distance, a polymeric mesh layer has been printed for previously mentioned reasons. The blue inset of figure 3d shows the aforementioned electrodes and the red inset shows that aluminum electrodes intimately wrap the 3D plastic object and polymeric mesh. Note that, blistered shapes shown in the blue inset are due to the Su8 mesh under aluminum film. Indeed, Su8-based ink suffers from the wellknown coffee ring effect. However, no mechanical failure was observed showing the robustness of the water transfer printing technology. Moreover, for external acquisition, an anisotropic conductive film (ACF 3M ECATT 9706) was used to connect the capacitive touchpad to a silver ribbon cable printed on Kapton. The capacitive touchpad has been transferred to the convex side of the curved object. As shown in figure 3e, the finger slides on the opposite side, faced to the electrodes, changing initial capacitance value. Arduino development board, connected to a computer via serial USB communication, has been used to measure the capacitance value between each electrode. Capacitances (Cn), are firstly measured using microcontroller analog outputs/inputs. Time constants (τn) of capacitors discharge τn = RCn are then estimated, where R is the pull-up internal resistor of the microcontroller.

As shown in figure 3f, when the finger position (Fpn) faces a capacitor (Cn), capacitance value decreases. Cn is then the most impacted (red bars) by Fpn enabling finger position sensing. One can note that neighboring capacitors are also, to a lesser extent, impacted by the finger (grey bars), depending on how it is placed. Moreover, a pong game using capacitive sensor as a remote has been demonstrated using "Processing" software (see figure 3g and Video S3 in supplementary file).

In conclusion, we have demonstrated an industrially viable technology to intimately wrap thin films around 3D objects. Water transfer printing technology also named hydrodipping has been adapted from graphic arts to conformal electronics. The adapted process allows the fabrication of electronics on 2D substrate using conventional microelectronics materials and methods.

Films can be transferred to many kinds of surfaces including plastic, metal, woods, all those which can be dipped into water. Thin films have been successfully patterned at lithographical accuracy or by inkjet printing over a large-area showing the versatility of the concept. Metallic and polymeric thin films have been successfully wrapped around random 3D shapes whose materials show a wide span of wetting properties. Moreover, water transfer printing enables highly bended thin films because no rigid substrate is required during transfer step. Importantly, this technology has clearly solved a challenging issue in the field of conformal electronics. Indeed, patterns have been wrapped around a 3D holed object. The bent touchpad and metal/insulator/metal capacitors, demonstrate potential applications of this technology. We demonstrated that daily-life hand-sized objects are compatible with the aforementioned technique. We assume that the presented technique is only limited by the materials depositions equipment scale in view of a larger area processing.

Figure 1 :

 1 Figure 1: Patterned thin films wrapping 3D object using water-transfer process.

Figure 2 :

 2 Figure 2: Process capabilities. Optical observations of patterns conformally transferred to: a)

Figure 3 :Figure S1 .

 3S1 Figure 3: Capacitive devices on curved objects. a) Optical picture of 3-layers MIM capacitors

Figure S3 .

 S3 Figure S3. Optical pictures showing silver ink printed on PVA substrate for large-area applications.

Figure S4 :

 S4 Figure S4: PVA dissolution behavior as function of process variants.

Figure S5 .

 S5 Figure S5. Electrodes transferred to curved glass.

Figure S6 :

 S6 Figure S6: Aluminum conformal wrapping around 50 cents euro coin.

Figure S7 :

 S7 Figure S7: Aluminum interconnects wrapped around human hair

ACKNOWLEDGMENTS

This work is supported by the European Union through the European Regional Development Fund (ERDF), and by the French region of Brittany (project: IMPRIM'). We thank, R. Beal (ISTIC Fablab) and S. Pernon (Université de Rennes 1 -ISTIC) for part dealing with 3D printing.

AUTHOR INFORMATION

Corresponding Author

* E-mail: maxime.harnois@univ-rennes1.fr

Notes

The authors declare no competing financial interest.

Author Contributions

The manuscript was written through contributions of all authors.

All authors have given approval to the final version of the manuscript.