Xavier Leroy
email: xavier.leroy@inria.fr

Michel Mauny
email: michel.mauny@inria.fr

Dynamics in ML

Objects with dynamic types allow the integration of operations that essentially require run-time type-checking into statically-typed languages. This paper presents two extensions of the ML language with dynamics, based on what has been done in the CAML implementation of ML, and discusses their usefulness. The main novelty of this work is the combination of dynamics with polymorphism.

Introduction

Static typing (compile-time enforcement of the typing rules for a programming language) is generally preferred over dynamic typing (the production of run-time tests to check the rules), since static typing reports type violations earlier, and allows for generation of more efficient code. However, one has to revert to dynamic typing for programs that cannot be recognized type-safe at compile-time. This situation often reveals weaknesses of the type system used. Dynamic typing could be avoided, then, by employing more advanced type systems. For instance, it had long been believed that generic functions (functions that can be applied to arguments of different types) can only be supported by dynamically-typed languages, such as Lisp, until the advent of polymorphic type disciplines, such as the one of ML, that permit static typing of such functions.

In contrast, there are programming situations that seem to require dynamic typing in an essential way. A first example is the eval function (and similar meta-level operations), that takes a character string, evaluates it as an expression of the language, and returns its value. The type of the returned value cannot be known at compile-time, since it depends on the expression given as argument. Another example is structured input/output. Some runtime systems provide an extern primitive, that takes an object of any type, and efficiently outputs a low-level representation of the object to persistent storage. The object can be read back later on, possibly in another process, by the intern primitive. The extern function can easily be typed in a polymorphic type system; but this is not the case for the intern function, since the type of its result depends on the contents of the file being read. In order to guarantee type safety, it is clear that the values returned by eval or by intern must carry some type information at run-time, and that this type information must be dynamically checked against the type expected by the context.

As demonstrated above, dynamic typing cannot be avoided for a few highly specific functions. But we would like to retain static typing for the huge majority of functions that can be typechecked at compile-time. What we need is a way to embed dynamic typechecking within a statically-typed language. The concept of objects with dynamic types (or dynamics, for short) is an elegant answer to this need. A dynamic is a pair of a value v and a type expression τ , such that value v has type τ . From the standpoint of static typing, all dynamics belong to the built-in type dyn. Type dyn represents those values that are self-described as far as types are concerned; that is, those values on which run-time type checking is possible.

Continuing the examples above, function eval naturally returns dynamics, so its static type is string → dyn. Similarly, intern has type io channel → dyn, and the extern function will be made to accept arguments of type dyn only, since the external representation of an object should now include its type.

Two constructs are provided to communicate between type dyn and the other types in the language. One construct creates dynamics by taking an object of any type and pairing it with its static type. The other construct checks the internal type of a dynamic against some static type τ , and, in case of success, gives access to the internal value of the dynamic with type τ .

In this paper, we consider the integration of dynamics, as described above, into the ML language [START_REF] Milner | The definition of Standard ML (revised)[END_REF]. The main novelty of this work is the combination of dynamics with a polymorphic type discipline. This combination raises interesting issues that have not been addressed yet. The main published references on dynamics have only considered first-order types [START_REF] Abadi | Dynamic typing in a statically-typed language[END_REF], or first-order types with subtyping [START_REF] Cardelli | Amber. In Combinators and Functional Programming Languages[END_REF]. Polymorphism gets mentioned in [START_REF] Abadi | Dynamic typing in a statically-typed language[END_REF], but very briefly and very informally. An early, unpublished work by Mycroft [13] is said to consider the extension of ML with dynamics, but we were unable to find this draft.

The two extensions of ML with dynamics we present here are not mere proposals. The simpler one has been fully integrated into the CAML system [START_REF] Weis | The CAML reference manual[END_REF], the ML implementation developed at INRIA, for more than two years. It has grown to the point of stability where dynamics are used inside the CAML system. The second, more ambitious extension was also extensively prototyped in CAML. This practical experience enables us to discuss the main implementation issues involved by dynamics. It also gives some hints on the practical usefulness of dynamics in an ML system, both for user-level programming and system-level programming.

The remainder of this paper is organized as follows. Section 2 presents a first extension of ML with dynamics. After an informal presentation, we formalize typing and evaluation rules for dynamics within a small subset of ML, and discuss type inference and compilation issues. Section 3 extends the system previously described with the ability to destructure dynamics (both the type part and the value part), and rebuild dynamics with the components of the structure. We adapt the typing and evaluation rules of section 2 to this extension. Section 4 discusses the practical usefulness of the two systems, based on some significant uses of dynamics in the CAML environment. Finally, we give a few concluding remarks in section 5.

This section describes dynamics as they are implemented in CAML release 2.6 and later [18, chapter 8].

Presentation

The new construct dynamic a is provided to create dynamics. This construct evaluates a, and pairs it with (the representation of) the principal type inferred for a. For instance, dynamic 1 evaluates to (1, int), and dynamic true to (true, bool). In any case, the expression dynamic a is of type dyn, without any mention of the internal type of the dynamic.

Objects with polymorphic types can be put in dynamics, provided their type is closed : none of the type variables free in their type should be free in the current typing environment. For instance, dynamic(function x → x) is perfectly legal, since the identity function has type α → α, and α is a fresh type variable, that does not appear anywhere else. The resulting dynamic value is (function x → x, ∀α. α → α). The explicit quantification over α emphasizes the fact that the internal type of the dynamic is closed, and suggests that the internal value is really polymorphic: it will be possible to use it with several types.

On the other hand, function x → dynamic x is rejected: dynamic x is typed in the environment x : α, where x does not have a closed type. In this case, it is problematic to determine at compile-time the exact type of the object put into the dynamic: static typing says it can be any instance of α, that is, any type. To correctly evaluate the function above, the actual type to which α is instantiated would have to be passed at run-time. Since polymorphic functions can be nested arbitrarily, this means that all polymorphic functions, even those that do not build dynamics directly, would have to take type expressions as extra parameters, and correctly propagate these types to the polymorphic functions they call. We would loose one major benefit of static typing: that run-time type information is not needed except inside dynamic objects. Such is the reason why dynamics are required to be created with closed types.

To do anything useful with a dynamic, we must gain access to its internal value, bringing it back to the statically-typed world. A run-time type check is needed at that point to guarantee type safety. This check must ensure that the internal type of the dynamic does match the type expected by the context. This operation is called coercion of a dynamic. Coercion is customarily presented as a special syntactic construct (the typecase construct in [START_REF] Abadi | Dynamic typing in a statically-typed language[END_REF]). This construct binds the internal value of the dynamic to some variable. It also handles the case where the run-time type check fails, and another coercion must be attempted, or an exception raised.

In ML, these two aspects, binding and failure handling, are already covered by the pattern-matching mechanism. Hence, instead of providing a separate coercion construct, we integrate dynamic coercion within pattern-matching. Namely, we introduce a new kind of pattern, the dynamic patterns, written dynamic(p : τ ). This pattern selects all dynamics whose internal value matches the pattern p, and whose internal type agrees with the type expression τ . For instance, here is a function that takes a dynamic, and attempts to print a textual representation for it:

let print = function dynamic(x : int) → print int x | dynamic(s : string) → print string s | dynamic((x, y) : int × int) →
print string "("; print int x; print string ", "; print int y; print string ")" | x → print string "?" For type matching, two behaviors can be considered. The first one is to require that the internal type of the dynamic is exactly the same as the expected type, up to a renaming of type variables. The other behavior is to also accept any dynamic whose internal type is more general than the expected type. For instance, dynamic [ ], whose internal type is ∀α. α list, matches the pattern dynamic(x : int list) with the latter behavior, but not with the former. We have retained the latter behavior, since it seems more coherent with the statically-typed part of the ML language (where e.g. the empty list can be used in any context that expects a list of integers).

Type patterns are allowed to require a polymorphic type, as in dynamic(f : α → α). This pattern matches any dynamic whose internal type is as general or more general than the type in the pattern (e.g. β → β, or β → γ). As a consequence of these semantics, identifier f can safely be used with any instance of the type α → α in the right-hand side of the pattern-matching:

function dynamic(f : α → α) → f f
The type matching semantics guarantee that f will be bound at run-time to a value that belongs to all type instances of the type scheme α → α. This is the only case in ML where a variable bound by a function construct can be used with several types inside the function body.

In the example above, the type variable α is not a regular pattern variable such as f: it is implicitly quantified universally by the dynamic pattern, and therefore cannot be instantiated during the matching process. For instance, the pattern dynamic(x : α list) only matches a dynamic of the polymorphic empty list, not any dynamic of any list. As a consequence, a type pattern τ more general than a type pattern τ will match less dynamics than τ , in contrast with regular ML patterns. This means that in a dynamic matching, the most general type patterns must come first. To catch polymorphic lists as well as integer lists, one must write In the latter definition, the second case would never be selected, since the first case also matches dynamics with internal type α list.

Syntax

We now formalize the ideas above, in the context of the core ML language, enriched with pattern-matching and dynamics. The syntax of the language is as follows:

τ ::= int [] α [] τ → τ [] τ × τ [] dyn p ::= x [] i [] (p, p ) [] dynamic(p : τ ) a ::= x [] i [] function p 1 → a 1 | . . . | p n → a n [] a a [] (a, a ) [] let p = a in a [] dynamic a
Type expressions (with typical elements τ , σ) are either atomic types such as the type int of integers; type variables α; function types; product types; or the type dyn of dynamics.

Patterns (typical elements p,q) are either variables x; integer constants i; pairs of patterns; or dynamic patterns dynamic(p : τ ).

Finally, for expressions (typical elements a,b), we have variables x; integer constants i; functions with pattern matching on the argument; function application; pair construction; the let binding with pattern matching on the first expression; and the dynamic construct to build dynamics.

Typechecking

The typing rules for this calculus are given in figure 1. Most of the rules are just those of the core ML language, revised to take pattern-matching into account in the function and let constructs. Two additional rules present the creation and coercion of dynamics.

The rules define the predicate E a : τ , meaning "expression a has type τ in the typing environment E". An auxiliary predicate, p : τ ⇒ E, is used, meaning "pattern p has type τ and enriches the type environment by E". Here, E stands for a sequence of typing assumptions of the format x : Σ, where Σ is a type scheme: a type expression τ with zero, one or more type variables α k universally quantified.

E ::= (x : ∀α 1 . . . α n .τ ) *
We write E(x) for the type associated to x in E. If x is bound several times, we consider the rightmost binding only. We write F V (τ ) for the free variables of type expression τ , and F V (E) for the union of the free variables of all type schemes in E. Finally, Clos(τ, V ) stands for the closure of type τ with respect to those type variables not in the set V . It is defined by:

Clos(τ, V ) = ∀α 1 . . . α n .τ where {α 1 , . . . , α n } is F V (τ ) \ V .
The Clos operator is extended pointwise to type environments.

The only new rules are rule 7, that deals with dynamic creation, and rule 11, that deals with dynamic coercion. Rule 7 says that the expression dynamic a has type dyn, provided that a has a closed type τ : none of the free variables of τ are free in the current typing environment E.

(1) E i : int . . τ n . Then, dynamic(p : τ ) binds the same variables to the same values. As described above, all type variables free in τ 1 . . . τ n can be generalized. Hence, we take that dynamic(p : τ ) binds x 1 . . . x n to values of types Clos(τ 1 , ∅) . . . Clos(τ n , ∅).

(2) E(x) = ∀α 1 . . . α n .τ E x : τ [α k ← τ k ] (3) p k : σ ⇒ E k E, E k a k : τ E function p 1 → a 1 | . . . | p n → a n : σ → τ (4) E a : σ → τ E b : σ E a b : τ (5) E a : σ E b : τ E (a, b) : σ × τ (6) E a : σ p : σ ⇒ E E, Clos(E , F V (E)) b : τ E let p = a in b : τ (7) E a : σ F V (σ) ∩ F V (E) = ∅ E dynamic a : dyn (8) x : τ ⇒ x : τ (9) i : int ⇒ (10) p : τ ⇒ E p : τ ⇒ E (p, p ) : τ × τ ⇒ E, E (11)
p : τ ⇒ E dynamic(p : τ ) : dyn ⇒ Clos(E, ∅)

Type inference

The type system presented above enjoys the principal type property, just as the one of ML. The principal type is computed by a trivial extension of the Damas-Milner type inference algorithm [START_REF] Damas | Principal type-schemes for functional programs[END_REF] to handle the dynamic construct and dynamic patterns. For the dynamic a construct, the only difficulty is to ensure that the type of a is closed. It would not be correct to infer the most general type τ for a, and fail immediately if some variables in τ are free in the current typing environment: these variables may become instantiated to monomorphic types later. Consider the expression (function x → dynamic x) 1. Assuming the function part of the application is typed before the argument, the function is given type α → dyn, and dynamic x appears to build a dynamic with non-closed type α. But when the application is typed, α gets instantiated to int, and we know that the dynamic is created with internal type int. Hence, the closedness check must be delayed to the end of type inference. The CAML implementation proceeds as follows: when typing dynamic a, all type variables that are free in the inferred type for a and in the current (12) e x ⇒ e(x) For dynamic patterns dynamic(p : τ ), the expected type τ is given explicitly in the pattern, so there is actually nothing to infer. We just check that the pattern p is of type τ , and record the (polymorphic) types of the variables bound by p. We have considered inferring τ from the pattern p and the right-hand side a of the pattern-matching, but this seems quite difficult, since variables bound by p can be used with several different types in a.

(13) e i ⇒ i (14) e function . . . | p k → a k | . . . ⇒ [e, . . . | p k → a k | . . .] (15) e a ⇒ [e , . . . | p k → a k | . . .] e b ⇒ v v < p k ⇒ e e , e a k ⇒ w k minimal e a b ⇒ w (16) e a ⇒ v e b ⇒ w e (a, b) ⇒ (v, w) (17) e a ⇒ v v < p ⇒ e e, e b ⇒ w e let p = a in b ⇒ w (18) e a ⇒ v T ype(a) = τ e dynamic a ⇒ dynamic(v : τ ) (19) v < x ⇒ x ← v (20) v < p ⇒ e v < p ⇒ e (v, v ) < (p, p ) ⇒ e, e (21) i
< i ⇒ (22) v < p ⇒ e θτ = σ dynamic(v : τ ) < dynamic(p : σ) ⇒ e

Evaluation

We now give call-by-value operational semantics for our calculus. Expressions are mapped to values, that is, terms with the following syntax: The evaluation rules are given in figure 2. They closely follow the structure of the typing rules. A first set of rules define the predicate e a ⇒ v, meaning "expression a evaluates to value v in environment e". The remaining rules define the auxiliary predicate v < p ⇒ e, meaning "value v matches pattern p, binding variables in p to values as described by e".

v ::= i [] (v, v ) [] dynamic(v : τ ) [] [e, p 1 → a 1 | . . . | p n → a n ] e ::= (x ← v) *
Since most rules are classical, we only detail the two rules dealing with dynamics. Rule 18 expresses that evaluating dynamic a amounts to evaluating a, and pairing its value with the static type of a. The type of a is not mentioned in the expression dynamic a, so there are some technicalities involved in defining precisely what it is. We assume all expressions a considered here are subterms of a given closed, welltyped term a 0 (the program). Let D be the principal type derivation for a 0 (the one that is built by the type inference algorithm). For each subterm a of a 0 , T ype(a) is defined as the type given to a in D. (In practice, dynamic expressions are annotated with their types during typing.)

Rule 22 defines the semantics of pattern matching over dynamics. The internal type τ of the dynamic is required to be more general than the type σ expected by the pattern: there must exist a substitution θ of types for type variables such that θτ is σ. The internal value of the dynamic is recursively matched against the value part of the dynamic pattern.

Compilation

In the current CAML implementation, internal types of dynamics are represented by the following term-like structure:

type gtype = Gvartype of int | Gconsttype of int × gtype list
Type constructors are identified by unique stamps instead of names to correctly handle type redefinition. Type variables are also encoded as integers. The code generated for dynamic a simply pairs the value of a with the structured constant representing T ype(a) as a gtype. For pattern matching, CAML provides a library function ge gtype, that takes two types and tests whether the first one is more general than the second one. The code generated for pattern matching on dynamics simply calls ge gtype with the internal type of the dynamic, and the expected type (again, a structured constant of type gtype). Only when ge gtype returns true is the sequence of tests matching the internal value against the pattern entered. Those tests were compiled assuming that the value being tested belongs to the expected type for the dynamic; therefore, it would be incorrect to match the internal value first, and then the internal type.

Little effort went into making run-time type tests faster. We have not yet encountered CAML programs that need to perform dynamic coercions inside tight loops. In case coercion speed becomes an important issue, we could first switch to the following representation for internal types of dynamics:

type gtype = Gvartype of gtype option ref | Gconsttype of int × gtype list
This representation makes it possible to perform instantiations by physical modifications on the type, which is more efficient than recording them separately as a substitution. (These physical modifications are undone at the end of the matching.)

Then, we could perform partial evaluation on the ge gtype predicate, since its second argument is always known at compile-time. Conventional pattern-matching compilation techniques [START_REF] Simon | The implementation of functional programming languages[END_REF] do not apply directly, however, since they consist in specializing term-matching predicates on their first argument (the more general term), not on the second one (the less general term). Specializing a matching predicate such as ge gtype on its second argument turns out to be just as hard as the more general problem of specializing a unification predicate on one of its arguments. The latter problem has been extensively studied in the context of Prolog compilation. A popular solution is the Warren Abstract Machine and its compilation scheme [START_REF] David | An abstract Prolog instruction set[END_REF][START_REF] Maier | Computing with logic: logic programming with Prolog[END_REF]. Most of the techniques developed there apply to our problem. We shall detail this issue at the end of section 3.6.

Non-closed types in dynamic patterns

This section presents an extension of the system presented above that makes it possible to match dynamic values against dynamic patterns with incomplete type information. This enables destructuring dynamics without specifying their exact type.

Presentation

With the previous system, the internal value of a dynamic can only be extracted with a fixed type. This turns out to be insufficient in some cases. Let us continue the print example of section 2.1. For product types, we would like to have a single case that matches all dynamics of pairs, prints the parentheses and comma, and recursively calls the print function to print the two components of the pair. This cannot be done with the system above: the pattern dynamic((x, y) : α × β) will only match dynamics whose internal type is at least as general as ∀α∀β. α × β, definitely not all dynamics whose internal type is a pair type. What we need is to have type variables in dynamic patterns that are not universally quantified, but rather existentially quantified, so that they can be bound to the corresponding parts of the internal type of the dynamic.

We now give a more complete version of the print function, with explicit universal and existential quantification for type variables in dynamic patterns. We will use it as a running example in this section.

Typing existential quantification

Let us detail first how these existentially quantified type variables behave when typing the right-hand side of the pattern-matching. Such a variable α can be bound to any actual type at run-time. Hence, at compile-time, we should make no assumptions about type α, and treat it as an abstract type. That is, type α does not match any type except itself; and type α must not escape the scope of the pattern-matching that binds it: α is not allowed to be free in the type of the returned value. As a consequence, the following two functions are rejected: function ∃α. dynamic(x : α) → x = 1 function ∃α. dynamic(x : α) → x while this one is perfectly legal:

function ∃α. dynamic((f, x) : (α → int) × α) → f x
and can be applied to dynamic(succ, 2) as well as to dynamic(int of string, "3").

There is one important difference between existentially bound type variables and abstract types: the actual type bound to such a type variable is available at run-time. Given an object a whose static type contains a variable α existentially bound, it is possible to build a dynamic from this object. The internal type of the dynamic will be the "true" type for a: its static type where the binding of α has been performed. Cases (3) and (5) in the print function illustrates this feature: when the matching with ∃α. dynamic(x :: l : α list) succeeds, two dynamics are created, dynamic x with internal type the type τ bound to α; and dynamic l with internal type τ list. This transforms a dynamic of a non-empty list into the dynamic of its head, and the dynamic of its tail, thus allowing recursion on the list.

Mixed quantifications

Existentially quantified variables can be freely mixed with universally quantified variables inside type patterns. Then, the semantics of the matching depends on the relative order in which these variables are quantified. This is illustrated by cases ( 7) and (8) in the print example -two modest attempts at printing functional values.

In case [START_REF] Gordon | Adding eval to ML[END_REF], the pattern is ∃α.∀β.dynamic(f : α → β). Since α is bound before β, variable α only matches type expressions that do not depend on β. For instance, a dynamic with internal type γ → γ is rejected. The functions selected by the pattern above are exactly those returning a value of type β for all β. Since no such value exists in ML, the selected functions never terminate (or always raise an exception), hence they are printed as function x → ⊥.

In case ( 8), the pattern is ∀α.∃β.dynamic (f : α → β). Here, β is bound after α; hence β can be instantiated by type expressions containing α. For instance, this pattern matches a dynamic with type γ → γ list, binding β to α list. This pattern catches a class of functions that operate uniformly on arguments of any type. These functions cannot test or destructure their argument, but only put it in data structures or in closures. Therefore, if we apply such a function to a symbolic name x, and recursively print the result, we get a representation of the function body, with x standing for the function parameter1 . (Actually, the printed function is extensionally equivalent to the original function, assuming there are no side-effects.)

In the presence of mixed quantification, the rules for typing the right-hand side of pattern-matchings outlined above have to be strengthened: it is not always correct to treat an existentially quantified type variable as a new abstract atomic type. Consider:

function ∀α.∃β.dynamic(f : α → β) → f 1 = f true
Assuming f : ∀α. α → β, the expression f 1 = f true typechecks, since both applications of f have type β. Yet, when applying the function above to dynamic (function x → x), the matching succeeds, f 1 evaluates to 1, f true evaluates to true, and we end up comparing 1 with true -a run-time type violation. Since the actual value of β is allowed to depend on α, static typechecking has to assume that β does depend over α, and treat two occurrences of β corresponding to different instantiations of α as incompatible. This is achieved by taking β in the right-hand side of the matching to be a type constructor parameterized by α. To avoid confusion, we shall write β for the type constructor associated to type variable β. Therefore, we now assume f : ∀α. α → β(α) for the typing of f 1 = f true, and this leads to a static type error, since the two sides of the equal sign have incompatible types β(int) and β(bool). However, f 1 = f 2 is well-typed, since both sides have type β(int). The general rule is: for the purpose of typing the right-hand side of a pattern-matching, existentially quantified type variables β are replaced by the type expression β(α 1 , . . . , α n ), where α 1 . . . α n is the list of those type variables that are universally quantified before β in the pattern. This transformation is known in logic as Skolemization.

Multiple dynamic matching

Type variables are quantified at the beginning of each case of the pattern-matching, not inside each dynamic pattern. This makes no difference for universally quantified variables. However, existentially quantified variables can be shared among several dynamic patterns, expressing sharing constraints between the internal types of several dynamics. For instance, the "dynamic function application" example of [START_REF] Abadi | Dynamic typing in a statically-typed language[END_REF] can be written as:

function ∃α.∃β. (dynamic(f : α → β), dynamic(x : α)) → dynamic(f x)
This function takes a pair of two dynamics, applies the first one (which should contain a function) to the second one, and returns the result as a dynamic. It ensures that the type of the argument is compatible with the domain type of the function.

Type variables can be shared among two dynamic patterns of the same matching; but we shall prohibit sharing between patterns belonging to different matchings (curried dynamic matching). In other terms, all cases in a pattern matching are required to be closed: all type variables contained in dynamic patterns should be quantified at the beginning of the corresponding matching. For instance, it is not possible to write the dynamic apply function in the following way (as it originally appears in [START_REF] Abadi | Dynamic typing in a statically-typed language[END_REF]):

function ∃α.∃β. dynamic(f : α → β) → function dynamic(x : α) → dynamic(f x)
This violates the requirement above, since α is bound by the outermost matching, and mentioned in the innermost one. The reasons for this restriction are mostly pragmatic: curried dynamic matching, in conjunction with polymorphic dynamics, can require postponing some type matching in an outer dynamic matching until an inner dynamic matching is performed. Our closedness condition on pattern matching cases rules out these nasty situations, without significantly reducing the expressive power of the language: curried dynamic application can still be written as

function df → function dx → match (df, dx) with . . .
at the expense of a later error detection, in case df is not a dynamic of function.

Syntax

The only syntactic change is the introduction of a sequence of quantifiers in front of each case in pattern matchings.

a ::= . . . [] function Q 1 p 1 → a 1 | . . . | Q n p n → a n Q ::= [] ∀α.Q [] ∃α.Q (23) Q k p k : σ ⇒ E k E, E k a k : τ F SC(τ ) ∩ BV (Q k ) = ∅ function Q 1 p 1 → a 1 | . . . | Q n p n → a n : σ → τ (24) Q x : τ ⇒ x : τ (25) Q i : int ⇒ (26) Q p : τ ⇒ E Q p : τ ⇒ E Q (p, p ) : τ × τ ⇒ E, E (27) F V (τ ) ⊆ BV (Q) Q p : τ ⇒ E θ = S( , Q) Q dynamic(p : τ ) : dyn ⇒ Clos(θτ, ∅)
Figure 3: Typing rules with explicit quantification in type patterns

We will always assume that variables are renamed so that quantifier prefixes Q never bind the same variable twice. We write BV (Q) for the set of variables bound by prefix Q.

Typechecking

We introduce the Skolem constants at the level of types. To each type variable α, we associate the type constructor α, with variable arity.

τ ::= . . . [] α(τ 1 , . . . , τ n ) Skolem constants are not permitted in the type part of dynamic patterns, nor in the internal types of dynamic values. We shall write τ 0 , σ 0 for type expressions free of Skolem constants. We define F SC(τ ), the free Skolem constants of type τ , as the set of all variables α such that type constructor α appears in τ .

The new typing rules for functions and for patterns are shown in figure 3. For each case Qp → a in a function definition, the pattern p is typed taking Q into account. The proposition p : σ ⇒ E now takes Q as an extra argument, becoming Q p : σ ⇒ E. The Q prefix is carried unchanged through all rules, and it is used only in the rule for dynamic patterns. There, in the types of all identifiers bound by the pattern, we replace existentially quantified type variables by the corresponding Skolem functions. This is performed by the substitution θ = S( , Q), defined inductively on Q as follows:

S(α 1 . . . α n , ) = id S(α 1 . . . α n , ∀α.Q) = S(α 1 . . . α n α, Q) S(α 1 . . . α n , ∃α.Q) = [α ← α(α 1 , . . . , α n )] • S(α 1 . . . α n , Q)
Typing of action a proceeds as previously. We simply check that the type of τ does not contain any Skolem constants corresponding to variables bound by Q.

(28) e int ⇒ int (29) e α ⇒ α (30) e dyn ⇒ dyn

(31) e σ ⇒ σ 0 e τ ⇒ τ 0 e σ → τ ⇒ σ 0 → τ 0 (32) e σ ⇒ σ 0 e τ ⇒ τ 0 e σ × τ ⇒ σ 0 × τ 0 (33) e(α) = [α 1 . . . α n ].τ 0 e τ k ⇒ τ 0 k e α(τ 1 . . . τ n ) ⇒ τ 0 [α k ← τ 0 k ] (34) e τ ⇒ τ 0 e a ⇒ v e dynamic a ⇒ dynamic(v : τ 0 ) (35) e a ⇒ [e 1 , . . . | Q k p k → a k | . . .] e b ⇒ v k minimal Q k v < p k ⇒ e 2 ; Γ Solve(Q k , Γ) = e 3 e 1 , e 2 , e 3 a k ⇒ w e a b ⇒ w (36) Q v < x ⇒ x ← v ; (38) Q v < p ⇒ e ; Γ Q v < p ⇒ e ; Γ Q (v, v ) < (p, p ) ⇒ e, e ; Γ, Γ (37)
Q i < i ⇒ ; (39) Q v < p ⇒ e ; Γ F V (τ ) ∩ BV (Q) = ∅ Q dynamic(v : τ ) < dynamic(p : σ) ⇒ e ; Γ, τ = σ

Evaluation

The introduction of existential type variables in dynamic patterns significantly complicates the semantics of the language, both for dynamic creation, and for dynamic matching. The modified evaluation rules are shown in figure 4.

For dynamic creation (rule 34), the evaluation of dynamic a now has to transform the static type τ inferred for a before pairing it with the value of a. Skolem constants representing existentially bound type variables are replaced by the actual types bound to these variables, properly instantiated. These bindings of type variables are recorded in the evaluation environment e. Hence the new syntax for evaluation environments:

e ::= (x ← v [] α ← [α 1 , . . . , α n ]τ 0 ) *
Since existential type variables may depend upon universal variables, existential variables are actually bound to a type context (a type expression with holes) instead of a simple type expression. We write type contexts as [α 1 , . . . , α n ]τ 0 , where type variables α 1 . . . α n are names for the holes. Rules 28-33 define the evaluation relation on types e τ ⇒ τ 0 , mapping a type expression τ to a type expression τ 0 without Skolem constants.

For dynamic matching during function application (rule 35), it is not possible anymore to perform dynamic type matching separately for each dynamic pattern, since patterns may share existentially quantified variables. Therefore, all dynamic type constraints are collected first, as a set of equations τ = σ, where τ is the internal type of a dynamic, and σ a type pattern. The pattern-matching predicate becomes Q v < p ⇒ e ; Γ, where Γ is the sequence of equations between types described above, and Q is the quantifier prefix for the matching (rules 36-39). The Q prefix is used in rule 39 to rename the internal types of dynamics, if necessary, so that their free type variables are not bound by Q. In a second phase, the function Solve is called to resolve the equations on types Γ, taking prefix Q into account. The precise definition of Solve is postponed to the next section. When the type matching succeeds, Solve returns the correct bindings for existentially quantified type variables. Then, evaluation of the right-hand side of the matching proceeds as usual.

Unification

The run-time matching between type patterns and internal types of dynamics amounts to a certain kind of unification problem, called unification under a prefix. This problem is studied extensively in [START_REF] Miller | Unification under a mixed prefix[END_REF], though in the very general setting of higher-order unification, while we only deal with first-order terms here. The first-order problem also appears in [START_REF] Miller | Lexical scoping as universal quantification[END_REF]. In our case, the problem consists in checking the validity of propositions of the format

q 1 α 1 . . . q n α n . σ = τ,
where the q k are either universal or existential quantifiers, and σ, τ are first-order terms of a free algebra. Unification under mixed prefix generalizes the well-known matching problem ("given two terms σ, τ , find a substitution θ such that θσ = τ ") and the unification problem ("given two terms σ, τ , find a substitution θ such that θσ = θτ "): writing α 1 . . . α n for the variables of σ, and β 1 . . . β n for the variables of τ , the matching problem is equivalent to

∀β 1 . . . ∀β n ∃α 1 . . . ∃α n . σ = τ,
and the unification problem to

∃β 1 . . . ∃β n ∃α 1 . . . ∃α n . σ = τ.
For the purpose of dynamic matching, we not only want to know whether the problem Q. σ = τ is satisfiable (Q is a quantifier prefix), but also to find minimal assignments for the variables existentially quantified in Q that satisfy the proposition. From now on, we shall treat variables universally bound in Q as constants. That is, we add such variables as term constructors with arity zero to the initial signature (int and dyn of arity zero, → and × of arity two).

Definition 1 A substitution θ is a Q-substitution iff for all variables α, all constants β contained in the term θα are bound before α in prefix Q. Definition 2 A substitution θ is a Q-unifier of σ and τ iff θσ = θτ , and θ is a Q-substitution. If such a substitution exists, σ and τ are said to be Q-unifiable.

Proposition 1 The proposition Q. σ = τ is satisfiable if and only if σ and τ are Q-unifiable.

Proposition 2 Two terms σ and τ are Q-unifiable if and only if σ and τ are unifiable, and their most general unifier is a Q-substitution.

Proof:

The "if" part is obvious. For the "only if" part, let θ be a Q-unifier of σ and τ . Since θσ = θτ , the terms σ and τ are unifiable. Let µ be their most general unifier. Let φ be a substitution such that θ = φ•µ. For all variables α, the constants contained in µα are a subset of those contained in θα. Since θ is a Q-substitution, all constants in µα are also bound before α in Q. Hence µ is a Q-substitution. 2

This result trivially gives an algorithm to compute the most general Q-unifier of σ and τ : compute the most general unifier of σ and τ , using Robinson's algorithm, and check that it is a Q-substitution.

We can now define the function Solve used in evaluation rule 32. It takes a prefix Q and a set Γ of equations σ 1 = τ 1 . . . σ n = τ n . Since Q binds the variables in the τ k only, prefix Q is first completed to bind the variables in the σ k also. Let α 1 . . . α m be the variables in the σ k . We take Q = Q.∃α 1 . . . ∃α m . (None of the α k is bound by Q, and they can be instantiated to any type.) Let µ be the most general Q -unifier of σ 1 ×. . .×σ n and τ 1 ×. . .×τ n , as computed by the algorithm above. Substitution µ is transformed into an evaluation environment, by adding bindings for the variables that are existentially quantified in Q. More precisely, we take Solve(Q, Γ) to be s(µ, , Q), where s is the run-time counterpart of the Skolemization function S used for static typing in section 3.3:

s(µ, α 1 . . . α n , ) = s(µ, α 1 . . . α n , ∀α.Q) = s(µ, α 1 . . . α n α, Q) s(µ, α 1 . . . α n , ∃α.Q) = α ← [α 1 , . . . , α n ]µα, s(µ, α 1 . . . α n , Q)

Compilation

The semantics given above are quite complicated, so it is no surprise their implementation turns out to be delicate. The main difficulty is unification under a prefix Q. Efficient algorithms are available for the regular unification phase. It remains to quickly check that the resulting substitution is a Q-substitution. This check can actually be integrated within the occur check, at little extra cost. The idea is to reflect dependencies by associating ranks (integers) to type variables. Variables bound by Q are statically given ranks 0, . . . , n from left to right. Other variables (i.e. those in the internal types of dynamics) are considered bound at the end of Q, and therefore given rank ∞. When identifying two variables α and β, the resulting variable is given rank min(rank(α), rank(β)). Then, binding existential variable α to a constructed type τ is legal iff:

1. (occur check) α does not occur in τ 2. (rank check) τ does not contain any universal type variable whose rank is greater than the rank of α.

As in the case of simple dynamics (section 2.5), the easiest way to implement type matching is to call at run-time a unification primitive, with the type pattern (annotated by rank information) as a constant argument. Partial evaluation of the unification primitive on the type pattern is desirable, not only to speed up type matching, but also to provide a cleaner handling of run-time type environments: after specialization, the bindings for the existential type variables could be recorded e.g. on the stack or in registers, as for regular variables; without specialization, the unification primitive would return a data structure containing these bindings, and less efficient code would be generated to access these bindings.

Specializing unification on one of its arguments is not much harder than specializing matching on its second argument (section 2.5). The techniques developed for the Warren Abstract Machine [START_REF] David | An abstract Prolog instruction set[END_REF][START_REF] Maier | Computing with logic: logic programming with Prolog[END_REF] directly apply, with the exception of the extra rank check. For instance, the WAM does not perform occur check for the initial binding of an existential variable, while we have to check ranks even in this case. Another difference is that backtracking is always "shallow", in the WAM terminology, since ML pattern-matching is deterministic. This simplifies the handling of the trail.

During the summer of 1988, the first author integrated a prototype unification compiler in the CAML system, following the ideas above. The CAML patternmatching compiler was modified to implement unification semantics as well as matching semantics, depending on flags put on the patterns. This low-level mechanism allowed performing unification on some parts of a data structure, and regular patternmatching on the other parts. Then, dynamic patterns dynamic(p : τ ) were simply expanded after type inference into product patterns (p, repr(τ )), where repr(τ ) is the pattern that matches all internal representations of types matching τ . The pattern repr(τ ) is marked to use unification semantics.

The only missing feature from what we have described above was rank check. At that time, we considered only dynamic patterns where all universal type variables come first, followed by all existential variables. Rank check could have been added with little modifications.

Dynamic matching benefited from all optimizations performed by the patternmatching compiler, including factorization of tests between cases, and utilization of typing informations. As a consequence, dynamic matching was performed quite efficiently. However, we agreed that this efficiency was not worth the extra complication of the compiler, and this prototype was never carried to the point it could be released.

Assessment

This section discusses the practical usefulness of the two propositions above, drawing from our practical experience with the CAML system.

Interfacing with system functions

Dynamics makes it possible to provide an interface with a number of system functions that cannot be given a static type in ML. Without dynamics, these functions could not be made available to the user in a type-safe way. In the CAML system, these functions include:

• extern : extern channel × dyn → unit and intern : intern channel → dyn, to efficiently write and read data structures on persistent storage, preserving sharing inside the structure. A typical use is, for a separate compiler, to communicate compiled object code with its linker, and to save and reload symbol tables representing compiled module interfaces.

• eval syntax : ML → dyn, to typecheck, compile, and evaluate a piece of abstract ML syntax (type ML). This makes it easy to provide CAML as an embedded language inside a program. For instance, the Calculus of Construction [START_REF] Huet | The Calculus of Constructions, documentation and user's guide[END_REF], a proof development environment, provides the ability to interactively define proof tactics written in CAML, and to apply them on the fly. The CAML macro facility [18, chapter 18] also makes use of eval syntax, since a macro body is an arbitrary CAML expression whose evaluation leads to the substituted text.

• MLquote : dyn → ML, which is one of the constructors of the datatype representing abstract syntax trees. This constructor embeds constants of arbitrary types inside syntax trees. These constants are produced by compile-time evaluation (e.g. macro expansion and constant folding).

• print : dyn → unit, to print a dynamic value in ML syntax. CAML cannot provide a polymorphic printing function with type α → unit, due to some optimizations in the data representation algorithm, that makes it impossible to decipher the representation of a data without knowing its type.

In these examples, the returned dynamics are generally coerced to fully known types, usually monomorphic. Therefore, we do not see the need for existential type variables there, and the simpler dynamic system presented in section 2 seems largely sufficient. In practice, the restriction encountered first is not that dynamics can only be coerced to closed types, but that dynamics can only be created with closed types. This prevents the print function from being called by a polymorphic function to print its polymorphic argument, for instance. This is often needed for debugging purposes.

Ad-hoc polymorphism

ML polymorphism is uniform: polymorphic functions operate in the same way on arguments of several types. In contrast, ad-hoc polymorphism consists in having generic functions that accept arguments of several types, but operate differently on objects of different types. Prime examples are the print function or the equal predicate: different algorithms are used to print or compare integers, strings, lists, or references. Several extensions of functional languages have been proposed, that support the definition of such generic functions, including type classes [START_REF] Wadler | How to make ad-hoc polymorphism less adhoc[END_REF] and runtime overloading [START_REF] Rouaix | Safe run-time overloading[END_REF].

Dynamics provide a naive, but easy to understand, way to define generic functions. As demonstrated above in the print example, dynamics permit joining predefined functions on atomic types (print int, print string) and functions on data structures (pairs, lists), that recurse on the components of the structures -the main operation in defining generic functions. Another important aspect of generic functions is extensibility: whenever a new datatype is defined, they can be extended to deal with objects of the new type as well. This can also be supported in the dynamic implementation, by keeping in a reference a list of, say, printing functions with type dyn → unit, to be applied until one succeeds whenever none of the standard cases apply. It should be pointed out that this implementation of generic functions with dynamics has several major drawbacks. First, because of the restrictions on dynamic creation, polymorphic functions that need to call print have to take dynamics themselves. This is not too serious for print, but would be prohibitive for heavily used generic functions, such as equal: all functions on sets, association lists, etc., would have to operate on dynamics, thus dramatically reducing accuracy of static typing and efficiency of compiled code. Moreover, nothing statically prevents the print function from being applied to objects that have no printing method defined. This important class of type errors will only be detected at run-time. Finally, such an implementation of generic functions is rather inefficient, since dynamics are built and coerced at each recursive call.

Type classes and run-time overloading techniques seem more realistic in this respect. They statically guarantee that generic functions can only be applied to objects on which they are defined. They perform type matching at compile-time whenever possible. And run-time type information can usually be arranged as dictionaries of methods, allowing faster method selection than dynamic type matching.

We have presented two extensions of ML with dynamic objects. The simpler one has proved quite successful for interfacing user code with some important system functions in a type-safe way. Its implementation cost remains moderate. The other extension, that generalizes the dynamic patterns to include both universal and existential variables in the type part, makes it possible to work on dynamics without coercing them to fixed types. Its semantics are more delicate, and therefore harder to implement. We lack strong evidence of its practical usefulness. We have presented one promising application: writing generic functions such as print in a way that is conceptually simpler than type classes. However, the usability of these functions is limited by the restriction that dynamics must be created with closed types. This restriction can be lifted, either by passing type information at run-time to all polymorphic functions, or by examining the call chain at dynamic creation time to reconstruct the instantiations of type variable -a technique developed for tagless garbage collection [START_REF] Appel | Run-time tags aren't necessary[END_REF][START_REF] Goldberg | Tag-free garbage collection for strongly typed programming languages[END_REF]. It remains to estimate the run-time penalty incurred.

  function dynamic(x : α list) → . . . | dynamic(x : int list) → . . . instead of the more intuitive definition function dynamic(x : int list) → . . . | dynamic(x : α list) → . . .

Figure 1 :

 1 Figure 1: Typing rules

Figure 2 :

 2 Figure 2: Evaluation rules

A

  value is either an integer i; a pair of values; a dynamic value dynamic(v : τ ) (a pair of a value v and a type expression τ ); or a closure [e, p 1 → a 1 | . . . | p n → a n ] of function body p 1 → a 1 | . . . | p n → a n by evaluation environment e. Evaluation environments map identifiers to values in the same way as typing environments map identifiers to types. For dynamic values, all type variables in the type part are considered universally quantified there; hence, two dynamic values are identified up to a renaming of their type variables.

  type fun arg = Arg of string in let rec print = function dynamic(i : int) → (1) print int i | dynamic(s : string) → (2) print string "\""; print string s; print string "\"" | ∃α.∃β.dynamic((x, y) : α × β) → (3) print string "("; print(dynamic x); print string ", "; print(dynamic y); print string ")" | ∃α.dynamic([ ] : α list) → (4) print string "[ ]" | ∃α.dynamic(x :: l : α list) → (5) print(dynamic x); print string " :: "; print(dynamic l) | ∀α.dynamic(f : α → α) → (6) print string "function x → x" | ∃α.∀β.dynamic(f : α → β) → (7) print string "function x → ⊥" | ∀α.∃β.dynamic(f : α → β) → (8) let s = gensym() in print string "function "; print string s; print string " → "; print(dynamic(f (Arg s))) | dynamic(Arg(s) : fun arg) → (9) print string s | ∃α.∃β.dynamic(f : α → β) → (10) print string "function x → . . . " | d → (11) print string "?"

Figure 4 :

 4 Figure 4: Evaluation rules with explicit quantification in type patterns

exception

  Cannot print; ; type fun arg = Arg of string in let printers = ref ([ ] : (dyn → unit) list) in let rec print = function . . . | d → let rec try print = function [ ] → print string "?" | f :: rest → try f d with Cannot print → try print rest in try print !printers and new printer f = printers := f :: !printers Assuming, for instance, that type foo = A | B of int × foo has been defined, we could add a printer for type foo as follows: new printer (function dynamic(A : foo) → print string "A" | dynamic(B(x, y) : foo) → print string "B("; print (dynamic x); print string ", "; print (dynamic y); print string ")" | x → raise Cannot print)

To avoid any confusion between the formal parameter and constants mentioned in the function body, formal parameters are represented by a local type fun arg = Arg of string. This ensures that the given function cannot create any terms of type fun arg, unless it is the print function itself. Fortunately, the self-application print(dynamic print) selects case[START_REF] Miller | Lexical scoping as universal quantification[END_REF] of the definition.


E Cueto
email: ecueto@unizar.es

M Laso

F Chinesta

Meshless stochastic simulation of micro-macro kinetic theory models

We present in this paper a numerical technique for the stochastic simulation of molecular models of visco-elastic fluids based on kinetic theory. The technique is based on the use of meshless methods and allows for an updated Lagrangian description of the conservation equations. It makes use of Natural Neighbor Galerkin schemes, that allow for a proper geometrical description of the domain as it evolves. The presented technique is especially well suited for the numerical simulation of free-surface flows. In this way, model molecules are associated with nodal positions such that they are advected with material velocities. Problems associated with lack of molecules in certain elements, for instance, as encountered in the basic implementation of CONNFFESSIT approaches, are thus avoided.

We present examples of validation and also performance tests of this technique applied to FENE and reptation (Doi-Edwards) models.

Introduction

Models arising from kinetic theory are nowadays a useful technique for the derivation of complex constitutive equations for the numerical simulation of visco-elastic flows. In fact, the majority of constitutive equations used in continuum numerical simulations of such fluids can be derived from kinetic theory [START_REF] Keunings | Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory[END_REF]. However, the derivation of a constitutive equation from kinetic theory models often involves closure approximations, that can have significant impact on the final predictions of the model thus constructed [START_REF] Chiba | On the fiber orientation in steady recirculating flows involving short fibers suspensions[END_REF][27] [START_REF] Pruliere | Recirculating flows involving short fiber suspensions: Numerical difficulties and efficient advanced micro-macro solvers[END_REF].

Micro-macro techniques aimed at resolving this problem. Although much more computationally expensive, these techniques, that couple coarse-grained molecular scale of kinetic theory to the macro scale of the continuum, have been used increasingly in the last years [START_REF] Laso | Calculation of viscoelastic flow using molecular models: the connffessit approach[END_REF] [START_REF] Gigras | Adaptive configuration fields: a new multiscale simulation technique for reptation-based models with a stochastic strain measure and local variations of life span distribution[END_REF]. In it, advantage is obtained by exploiting the equivalence of the Fokker-Planck equation, that describes the evolution of the probability distribution function governing the configuration state at the micro level, and an Itô stochastic differential equation. Instead of solving the deterministic Fokker-Planck equation, an equivalent stochastic differential equation is solved by means of a large enough ensemble of realizations of the stochastic process. The viscoelastic stress is thus obtained as an ensemble average of some function of the micro state at each molecule.

In the original CONNFFESSIT approach, the simulation begins by solving conservation equations, generally by using standard Eulerian or ALE Finite Element Methods. With the updated velocity field, a set of molecules are distributed on the flow, and their trajectory is computed in a Lagrangian framework. The Itô's stochastic differential equation is then integrated along each particle's trajectory and finally the viscoelastic stress is obtained by averaging the configuration state of the particles found at each finite element.

This approach presents, despite its success for many applications, some drawbacks (see the review by Keunings [START_REF] Keunings | Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory[END_REF] and references therein). First, molecules must be located within the finite element mesh. This tracking may be expensive for complex geometries. Second, a large number of particles must be located within each finite element in order to provide an accurate enough visco-elastic stress prediction, free of statistical noise.

These problems can be avoided if a Lagrangian strategy is implemented. In it, molecules are attached to the nodes, travelling together with material velocity. In this way, no element runs out of particles and no tracking algorithm is necessary to determine the element a particle is located in. The obvious problem associated with Lagrangian strategies is mesh distortion, or remeshing, if performed, and its associated numerical diffusion. However, in the middle nineties, the family of meshless methods [START_REF] Belytschko | Element-Free Galerkin Methods[END_REF] gained a tremendous popularity due precisely to their lack of sensitivity to mesh distortion. No mesh (in a FE sense) must be generated in meshless methods, since it is done by the method in a process transparent to the user.

Meshless methods thus opened the possibility of Lagrangian approaches in many branches of Computational Mechanics [START_REF] Martínez | Fixed mesh and meshfree techniques in the numerical simulation of injection processes involving short fiber suspensions[END_REF] [31] [1] [START_REF] González | A Natural Element updated Lagrangian strategy for free-surface Fluid Dynamics[END_REF]. In this paper we analyze how these methods can help to solve some of the problems associated to stochastic approaches to micro-macro models of visco-elastic fluid flow.

To this end, we have employed the Natural Element Method (NEM) [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]. This method offers some advantages over other meshless methods for this kind of applications. First of all, the enforcement of essential boundary conditions is straightforward if we employ the α-NEM [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] or C-NEM [START_REF] Yvonnet | A new extension of the Natural Element method for non-convex and discontinuous problems: the Constrained Natural Element method[END_REF] versions of the method. Second, tracking of free surface is also straightforward. α-NEM techniques notably simplify this issues such that no boundary marker or similar is necessary. Thus, the use of techniques such as Volume of Fluid (VoF) or similar, employed recently in the context of visco-elastic free surface flows [START_REF] Bonito | Numerical simulation of 3d viscoelastic flows with free surfaces[END_REF] can also be overcome.

The outline of the paper is as follows. In section 2 we describe the basics of the NEM, together with its most important properties. In section 3 we describe the models implemented, namely FENE and reptation models, and finally in section 4 we present some numerical results that demonstrate the potential of the method. The paper is completed with the conclusions.

The Natural Element Method

As mentioned before, one of the main advantages of meshless methods in general, and NEM in particular, is that no lack of accuracy is due to mesh distortion [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]. This is in sharp contrast with the Finite Element Method, whose behavior under distortion is well known many years ago [START_REF] Babuška | On the angle condition in the finite element method[END_REF].

In essence, the NEM is a Galerkin method, just like the FEM, but instead of piece-wise polynomial approximation, natural neighbor (and hence its name) approximation is used. In this section the construction of natural neighbor interpolation is reviewed, together with its most important properties. Prior to this definition we introduce some basic geometrical entities that are needed for further developments. each point within these regions is closer to the node to which the region is associated than to any other in the cloud. This kind of space decomposition is called a Voronoi diagram (also Dirichlet tessellation) of the cloud of points and each Voronoi cell is formally defined as (see figure 1):

Voronoi/Dirichlet diagrams

T I = {x ∈ R d : d(x, x I ) < d(x, x J ) ∀ J = I}, (1)
where d(•, •) is the Euclidean distance function.

The dual structure of the Voronoi diagram is the Delaunay triangulation1 , obtained by connecting nodes that share a common (d-1)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there being some so-called degenerate cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in 2D, as depicted in Fig. 1 (right)). Another way to define the Delaunay triangulation of a set of nodes is by invoking the empty circumcircle property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called natural neighbours and hence the name of the technique.

Equivalently, the second-order Voronoi diagram of the cloud is defined as

T IJ = {x ∈ R d : d(x, x I ) < d(x, x J ) < d(x, x K ) ∀ J = I = K}. (2)
Based on these definitions, different natural neighbour interpolation schemes have been proposed. We focus on Thiessen and Sibson interpolation [START_REF] Sibson | A Vector Identity for the Dirichlet Tesselation[END_REF]. The interested reader can look for other related interpolation schemes in [START_REF] Sukumar | Natural Neighbor Galerkin Methods[END_REF] and references therein.

Thiessen interpolation

The simplest of the natural neighbour-based interpolants is the so-called Thiessen's interpolant [START_REF] Thiessen | Precipitation averages for large areas[END_REF]. Its interpolating functions are defined as

ψ I (x) = 1 if x ∈ T I 0 elsewhere (3)
The Thiessen interpolant is a piece-wise constant function, defined over each Voronoi cell. It defines a method of interpolation often referred to as nearest neighbour interpolation, since a point is given a value defined by its nearest neighbour. Although it is obviously not valid for the solution of second-order partial differential equations, it can be used to interpolate the pressure in formulations arising from Hellinger-Reissner-like mixed variational principles, as proven in [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF].

Sibson's interpolation

The most extended natural neighbour interpolation method, however, is the Sibson interpolant [36] [37]. Consider the introduction of the point x in the cloud of nodes. Due to this introduction, the Voronoi diagram will be altered, affecting the Voronoi cells of the natural neighbours of x. Sibson [START_REF] Sibson | A Vector Identity for the Dirichlet Tesselation[END_REF] defined the natural neighbour coordinates of a point x with respect to one of its neighbours I as the ratio of the cell T I that is transferred to T x when adding x to the initial cloud of points to the total volume of T x . In other words, if κ(x) and κ I (x) are the Lebesgue measures of T x and T xI respectively, the natural neighbour coordinates of x with respect to the node I is defined as

φ I (x) = κ I (x) κ(x) . (4)
In Fig. 2 the shape function associated to node 1 at point x may be expressed as

φ 1 (x) = A abf e A abcd . (5)

Properties of Sibson's interpolation

Sibson interpolants have some remarkable properties that help to construct the trial and test functional spaces of the Galerkin method (see [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF] for proofs of the following properties).

x 1 Besides properties like continuity and smoothness (everywhere except at the nodes), Sibson interpolants possesses linear completeness (i.e., exact reproduction of a linear field).

Sibson interpolants can also reproduce linear functions exactly along convex boundaries. This is in sharp contrast to the vast majority of meshless methods. In addition, in [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] [14] [START_REF] Yvonnet | A new extension of the Natural Element method for non-convex and discontinuous problems: the Constrained Natural Element method[END_REF] distinct methods of imposing linear displacement fields along non-convex boundaries were developed. These are based on the use of α-shapes, ε-samplings or visibility criteria, respectively. So, essential boundary conditions can be imposed directly, as in traditional Finite Element methods.

In [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] and [START_REF] Cueto | Modeling three-dimensional piecewise homogeneous domains using the α-shape based Natural Element Method[END_REF] it was demonstrated that the construction of the Sibson interpolant over an α-shape [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] of the domain allows us to accurately extract the shape of the domain, defined in terms of nodes only, while ensuring linear interpolation along any kind of boundaries (convex or not). This property was later generalized for arbitrary clouds of points and a explicit definition of the domain through CAD techniques in [START_REF] Cueto | On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods[END_REF].

In the next section we study the implication of α-shapes in the development of the method here proposed. mensions this technique becomes prohibitive, since the markers should be boundary triangles.

In meshless methods models are constructed by a set of nodes only, and thus boundary tracking is usually performed by employing different strategies, in absence of such information. In particular, we have employed shape constructors to perform this task [START_REF] González | A Natural Element updated Lagrangian strategy for free-surface Fluid Dynamics[END_REF]. Shape constructors are geometrical entities that transform finite point sets into a multiply connected shape in general. Due to their importance in many areas, they have attracted much attention in Computational Geometry in the last years. In particular, we employ α-shapes [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF]. α-shapes define a one-parameter family of shapes S α (being α the parameter), ranging from the "coarsest" to the "finest" level of detail. α can be seen, precisely, as a measure of this level of detail. This means that all details of size less than α will be ignored in the geometry of the domain. This is a common practice in many conversions from CAD models (very often geometrically "exact") to FEM models, in which very small features of the geometry are ignored from a mechanical point of view.

An α-shape is a polytope that is not necessarily convex nor connected, being triangulated by a subset of the Delaunay triangulation of the points. The Delaunay triangulation is chosen as a departure point due to its uniqueness (other than degenerate cases). But the Delaunay triangulation is always defined on the convex hull of the cloud of nodes, and thus no associated shape is obtained. A filtering is then applied to the list of triangles (tetrahedra) on the triangulation. They are ordered according to the radius of their circumcircle (circumsphere).

If a level of detail α is then fixed, it is clear that, if the nodal sampling is dense enough -for this, taking h, the nodal spacing parameter, lower than α will suffice-, all triangles whose cirumcircle is bigger than α must be eliminated from the true geometry of the domain. Crudely speaking, they hide concave features of the geometry, by assumption bigger than α.

The theory of α-shapes was soundly established in [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF] for points in the plane and later in [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] for three-dimensional clouds. α-shapes constitute a provable geometry reconstruction algorithm, in the sense that there exist formal geometrical and topological guarantees on the final result, based on the size of the sampling, as mentioned above. The interested reader should consult [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] for more details.

In Fig. 3 an example of the just introduced theory is presented. It represents some instances of the finite set of shapes for a cloud in a intermediate step of the simulation of a wave breaking in a beach. Note how as α increases from Fig. 3(a) to (f), we obtain different shapes ranging from the cloud of nodes itself (for α = 0) to the convex hull of the cloud. For initial value problems, in which we know the starting configuration of the cloud, and thus the exact geometry of the domain at the beginning of the simulation, but not the subsequent geometries being described for different time steps, we employ a related strategy. We first fix the nodal sampling (so to speak, the nodal spacing parameter h, according to the needs of the simulation) and then fix α > h. All details of size slightly bigger, but close to, h will then be reproduced, as in Fig. 3(d). More details on this technique can be found in [START_REF] Galavís | Domain tracking in meshless simulations of free-surface flows[END_REF] and [START_REF] González | A Natural Element updated Lagrangian strategy for free-surface Fluid Dynamics[END_REF]. Obviously, in large geometrical transformations an α-adaptivity can be envisaged [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF].

Mechanical modelling

Following Keunings [START_REF] Keunings | Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory[END_REF], kinetic theory possesses three main ingredients. One is the configuration space, denoted here by X. It represents the set of all possible values for molecular configurations. The second ingredient characterizes the state of all the molecules at a given point in the model, through a probability distribution function, ψ(X, x, t). Under equilibrium conditions we assume that this distribution function takes the form ψ eq . In these conditions, Brownian and elastic forces are, of course, in equilibrium. Non-equilibrium situations provoke in the fluid the appearance of a viscoelastic stress, denoted by τ . Velocity and stress fields are linked by the following conservation equations, in absence of inertia terms:

Divσ = 0 (6) Divv = 0 (7) σ = -pI + τ + 2η γ ( 8
)
where γ represents the strain rate tensor, γ = 1 2 (∇v + ∇v T ), η the solvent viscosity (when it exists) and τ the polymer stress coming from the microscopic model. The resulting model defines a two-scales non-linear transient model because the non-linear and transient behavior of the microscopic description that we introduce later.

The simplest solution algorithm lies in the use of an explicit strategy that allows linearizing the microscopic model by decoupling both description scales. Thus, the weak formulation related to Eqs. ( 6) and ( 8) when the extra-stress τ is assumed known from the previous time step (τ n-1 ) writes:

Ω 2η γn : γ * dΩ - Ω p n I : γ * dΩ + Ω τ n-1 : γ * dΩ = 0 (9)
In absence of a solvent contribution, i.e. when η = 0 in Eq. ( 9), as encountered in usual models of entangled polymers and melts, a DEVSS formulations [START_REF] Baaijens | An Iterative Solver for the DEVSS/DG Method with Application to Smooth and Non-Smooth Flows of the Upper Convected Maxwell Fluid[END_REF] is retained.

To complete the formulation, kinetic theory provides also the equation governing the evolution of the probability distribution function ψ. This equation is known as the Fokker-Planck equation:

Dψ Dt = - ∂ ∂X • {Aψ} + 1 2 ∂ ∂X ∂ ∂X : {Dψ} (10)
where A represents a vector describing the drift exerted by the fluid on the function ψ, and D is a symmetric, positive-definite matrix accounting for brownian effects in the model. D/Dt represent material derivative. The expression, finally, that relates the obtained configuration state with its enforced state of stress is known as the Kramers formula:

τ = g(X)ψdX = g(X) (11)
where the brackets • denote an ensemble average over all the molecular conformational space at a physical point, and g is some function of the configuration state, depending on the particular model considered. Despite the existence of some incipient deterministic schemes able to solve efficiently Eq. ( 10) (see [START_REF] Ammar | On the Reduction of Kinetic Theory Models Related to Finitely Extensible Dumbbells[END_REF] [2] [3] [START_REF] Ryckelynck | On the a priori Model Reduction: Overview and recent developments[END_REF] [12] [START_REF] Poitou | On the solution of fokker-planck equations in steady recirculating flows involving short fiber suspensions[END_REF] and the references therein), stochastic simulation strategies are nowadays much more used, fact that justifies the choice retained in this work.

The stochastic approach to the problems takes benefit from the equivalence of the Fokker-Planck equation [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] to the following Itô's stochastic differential equation [START_REF] Öttinger | Stochastic Processes in Polymeric Fluids[END_REF]:

dX = Adt + B • dW (12)
where

D = B • B T (13)
and W represents a Wiener process. Equation [START_REF] Chinesta | On the reduction of stochastic kinetic theory models of complex fluids[END_REF] applies along individual molecule trajectories. Thus, in our approach, these molecules will be attached to the nodes, allowing for the integration of Eq. ( 12) by the method of characteristics, along the nodal paths. The integration is performed by means of the so-called Euler-Maruyama scheme:

X j n+1 = X j n + A(X j n , t n )∆t + B(X j n , t n ) • ∆W j n ( 14
)
where n refers to the current time step and j to the individual molecule being integrated.

In this paper two different models based on the stochastic approach to the Fokker-Planck equations have been implemented. The first one is the well-known FENE model for polymer solutions and the second one is the Doi-Edwards model of entangled polymers based on the repatation concept.

FENE model

The Finite Extension Non-linear Elastic (FENE) model of polymers represents a set of elastic dumbbells that do not interact with each other, flowing in a newtonian solvent. Each dumbbell is composed by two Brownian beads joined by a spring. Thus, the conformation space is given by X = Q, Q being the vector joining the beads. This spring is assumed to have a nonlinear behavior, such that it has a finite, limited extension Q 0 . Thus, the constitutive equation for these springs will be [START_REF] Owens | Computational Rheology[END_REF]:

F c = H 1 -Q 2 /Q 2 0 Q (15)
with H the spring constant. Each bead is subjected to an hydrodynamic drag given by

F h = -ζ [ ṙ -(κ • r + v 0 )] ( 16
)
where r represents the position vector of each bead, κ the gradient of the velocity field and then κ • r + v 0 denotes the undisturbed fluid velocity at the bead position (homogeneous deformations are assumed at the molecule scale). Finally, ζ represents the viscous friction coefficient existing between the beads and the solvent. Moreover, the beads are subjected to Brownian forces originated by the bombardment of solvent molecules

F b = -kT ∂ ln Ψ ∂Q (17)
The dimensionless stress tensor for the FENE model is given by

τ = -nk B T Q Q 1 - Q2 b + I (18)
where, in turn, Q = Q/(k B T /H) 1/2 , k B represents the Boltzmann constant. b represents the so-called "dimensionless finite extensibility parameter" given by b = HQ 2 0 /k B T . The dimensionless Fokker-Planck equation for this particular model reads [START_REF] Herrchen | A detailed comparison of various FENE models[END_REF]:

∂ ψ ∂ t = - ∂ ∂ Q • κ • Q ψ - 1 2 1 1 - Q2 /b Q ψ + 1 2 ∂ ∂ Q • ∂ ∂ Q ψ ( 19
)
In this last equation, κ represents the dimensionless form of the velocity gradient and t = t/λ H with λ H = ζ/4H a characteristic time of the polymer.

The dimensionless distribution function ψ is given by

ψ( Q, t) = ψ(Q, t) k B T H 3 2 (20)
The equivalent Itô's stochastic differential equation for this model, Eq. ( 12), will be given by

d Qt = κ • Qt - 1 2 1 1 - Q2 t /b Qt d t + dW t (21)
In this case, W t is a three-dimensional Wiener process. To numerically integrate this equation we employed the second-order semi-implicit predictorcorrector algorithm proposed in [START_REF] Öttinger | Stochastic Processes in Polymeric Fluids[END_REF] and also employed in [START_REF] Herrchen | A detailed comparison of various FENE models[END_REF]. Further details can be found in these references.

Doi-Edwards model

The reptation model implemented in this work is essentially the ubiquitous Doi-Edwards or Curtiss-Bird model, see [START_REF] Öttinger | Stochastic Processes in Polymeric Fluids[END_REF] and references therein. In it, the configuration space X is defined by a tube segment orientation vector u and the normalized contour label s ∈ [0, 1]. Thus, the probability distribution function ϕ will now be function of u, s, x and t. The expression ϕ(u, s, x, t)duds gives the probability of finding at time t and position x a tube segment oriented between u and u + du and containing the chain segment [s, s + ds]. The Fokker-Planck equation associated with this model reads

Dψ Dt = - ∂ ∂u • [(I -uu) • κ • uψ] + 1 π 2 τ d ∂ 2 ψ ∂s 2 (22)
with τ d the so-called disengagement or reptation time [START_REF] Keunings | Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory[END_REF]. The constitutive Kramers equation has been taken in this work as

τ = G uu (23)
The Fokker-Planck equation is now equivalent to two different Itô stochastic differential equations. The first one, which is actually deterministic, reads

du = (I -uu) • κ • udt (24)
stating that the evolution of the tube orientation is governed by the drift of the surrounding fluid, through its velocity gradient, κ. The second one,

ds = 2 π 2 τ d dW (25)
governs the evolution of the chain segment escaping the tube, and is purely diffusive. Again, W represents a Wiener process.

5 Model integration

Functional approximation

The variational form of the macroscopic description of the flow kinematics, given by a Hellinger-Reissner-like principle, see Eq. ( 9), is approximated by choosing suitable interpolation schemes for both velocities and pressure fields. It is well-known that not all approximations give rise to a stable method. The well-known LBB condition [START_REF] Babuška | The Finite Element Method with Lagrange multipliers[END_REF] is recognized as the condition to be fulfilled by such an approximation.

In this work we have employed a mixed Sibson-Thiessen approximation for velocities and pressures, respectively. Although it is known that this kind of approximation does not verify the LBB condition for some, very specific, cases (see [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF]), it often renders stable results and is very convenient from the computational point of view. Other approximations are also possible within the NEM. For instance, in [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF] an enriched approximation that verifies the LBB condition was presented, although it includes one degree of freedom more for each node. Recently, natural neighbor approximations with quadratic or higher order consistency have been also presented, see [START_REF] González | Higher-order Natural Element Methods: towards an isogeometric meshless method[END_REF] [START_REF] Yvonnet | Natural element approximations involving bubbles for treating mechanical models in incompressible media[END_REF]. Both possibilities give a stable approximation for this kind of variational formulations.

Model discretization

Once an appropriate interpolation scheme has been chosen, the variational form of the problem given by Eqs. ( 6)-( 8) is approximated in a Bubnov-Galerkin framework,

v h (x) = n I=1 φ I (x)v I ( 26
)
p h (x) = n I=1 ψ I (x)p I . ( 27
)
where n represents the number of nodes considered in the approximation (natural neighbors of point x). The functions ψ I (x) and φ I (x) in this work represent some form of natural neighbour interpolation, as mentioned before (Thiessen and Sibson, respectively, for the choice made in this work). This leads to the following system of algebraic equations:

K G G T M v p = f 0 ( 28
)
where

K IJ = Ω B T I CB J dΩ (29)
G IJ = - Ω BT I ψ J dΩ (30)
M IJ = - 1 λ Ω ψ I ψ J dΩ (31)
f I = Γt φ I tdΓ ( 32
)
where t represents the applied tractions on the boundary and

BI = [φ I,1 (x) φ I,2 (x)] ( 33
)
B I =   φ I,1 (x) 0 0 φ I,2 (x) φ I,2 (x) φ I,1 (x)   ( 34
)
C IJKL = η(δ IK δ JL + δ IL δ JK ) (35)
In the Eq. ( 33), φ I,1 and φ I,2 represent the derivatives of the shape function with respect to coordinates 1 and 2, i.e., x and y, respectively. It must be noted that, if we consider totally incompressible situations, M = 0.

An overview of the solution algorithm

The proposed procedure implies some modifications to the standard CON-NFFESSIT technique. The standard algorithm, as proposed by [START_REF] Laso | Calculation of viscoelastic flow using molecular models: the connffessit approach[END_REF] was composed of the following steps:

• To this end, a large enough number of particles are distributed on the geometry of the flow and convected with it. Keunings [START_REF] Keunings | Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory[END_REF] cites some difficulties associated to this procedure, such as the tracking of particles. Simple schemes for tracking particles can be highly inaccurate. In addition, it can be highly CPU demanding, due to the large number of particles needed to avoid numerical noise (in the order of 10 5 -10 7 ). Finally, to avoid statistical problems, each finite element must contain a minimum number of particles (10 2 -10 3 ), and this could not happen, depending on the conditions of the flow.

The proposed algorithm, however, proceeds as follows:

• Step 1: Solve the macro equations via NEM and obtain the velocity field v

• Step 2: Particles' path is nodal path, so advancing their position is straightforward The advantage of this procedure is two-fold. First, no tracking algorithm is necessary, since particles are attached to the nodes. Second, no element runs out of particles. The price to pay is to obtain the velocity gradient at nodal position, in order to perform Step 3 below. This is addressed in the next section.

•

Numerical issues

The correct implementation of the models described above requires to take into consideration the best way to compute numerically the velocity gradient, κ, at nodal positions. This is needed, since we assume that model's particles will be attached to the nodes and travel with material velocity along the characteristic lines. Three main possibilities arise, in view of the particular features of Sibson's interpolation, presented in section 2. This kind of interpolation is not differentiable at nodal positions (similarly to finite elements), and thus the velocity gradient is not defined at these positions.

One first possibility is to employ standard Gauss quadratures on Delaunay triangles and to project an average of these values to the nodes. To this end we employ the superconvergent patch recovery approach [START_REF] Zienkiewicz | The Superconvergent Patch Recovery and a posteriori error estimates. Part I: the recovery technique[END_REF], slightly modified to fit into NE methods. Let κ * be a matrix containing the desired nodal values of the velocity gradient at nodes neighboring a particular Delaunay triangle. We obtain a smoothed gradient field by interpolating these values through

κ * = Φκ * (36)
where Φ represents the matrix containing the NEM's shape functions values. The problem is now formulated as the minimization of the functional

Ω (κ * -κ h ) 2 dΩ (37)
with κ h the gradient field obtained by direct differentiation of the velocity field, and thus discontinuous. Here, by Ω, we denote the Delaunay triangle. After minimization,

κ * = Ω Φκ h dΩ Ω Φ T ΦdΩ (38)
We employ three-point Gauss quadrature on triangles. Since the number of natural neighbors of any point is usually more than three, the resulting system of equations is often ill-behaved. However, by employing conjugate gradient methods, for instance, in solving Eq. ( 38) we have found no problem, if the algorithm is feeded with an estimate of the result (the values at Gauss points, for instance, which should be close to the sought nodal values).

Once equation system [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF] has been solved, different values of κ * are predicted for each node belonging to neighboring triangles. The process finishes by averaging these values at each node.

This technique is perhaps the crudest one, but it is computationally less expensive than the other two, namely, the employ of stabilized nodal quadratures [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] [21] -which assumes a homogenized gradient value on each Voronoi cell-or the use of natural neighbor interpolants of continuity C 1 or higher. In this last case, the derivative of the velocity field would be obtained by simply differentiating the method's shape functions at the nodal positions. The interested reader can consult [START_REF] González | Higher-order Natural Element Methods: towards an isogeometric meshless method[END_REF] for more details in this last method.

We employed the first approach in order to minimize the computational cost of the global method, which is high by nature due to the stochastic character of the equations.

6 Numerical results

Validation of the proposed technique: start-up of shear flow

In order to validate the method proposed before, we have reproduced wellknown results from existing techniques. To this end, a homogeneous shear flow was considered given by v x = γy, v y = v z = 0, imposed at the nodes on the boundary of a squared domain.

The results obtained in [START_REF] Herrchen | A detailed comparison of various FENE models[END_REF] by employing a standard CONNFFESSIT approach for the FENE model were compared to ours. A cloud of 525 nodes was employed. A uniform velocity gradient was imposed at the integration points and the resulting stress field was analysed. The projection technique explained in Section 5.4 was employed. THe objective is thus to validate the projection algorithm, as well as the the numerical integration for the SDE.

Results for the evolution of the viscosity, first and second normal stress coefficients are shown in Fig. 4 for one particular node in the cloud. Results are similar for all nodes. In all cases, λ H • γ = 100.0 and b = 20.0.

In general, good agreement is observed with respect to the results report in the detailed study by Herrchen and Ötinger [START_REF] Herrchen | A detailed comparison of various FENE models[END_REF]. For long times, a somewhat more pronounced statistical noise is observed, due to fact that a limited number of stochastic realizations has been accomplished at each node. In this case, 10000 realizations were considered. Of course, this is not the best approach to simulate a homogeneous, rheometrical flow, but it was considered here as a means to validate the approach. More complex, non-homogeneous flows, with free surfaces, are considered in the following sections.

Die exit flow of a FENE fluid

In this section we consider the time evolution of a die swelling flow of a FENE fluid. The initial cavity and cloud of points is depicted in Fig. 5.

The flow of the fluid is forced by adding new nodes entering the channel with prescribed parabolic velocity. The maximum value of this inflow velocity was set to 0.1m/s and λ H = 1.0. These nodes incorporate particles with equilibrium probability distribution function ψ eq . Thus, the number of nodes in the model is increasing with time. The time step during the simulation, constant, was set to 0.005s, both for the flow kinematics solver and the Euler-Maruyama integration scheme for the SDEs. This was so for simplicity, but different strategies can also be implemented. For instance, updating nodal positions only after some time steps in the integration of the SDEs is also possible, and sometimes even convenient, given the small time increment needed for a proper result in the Euler-Maruyama scheme.

In Fig. 6 a snapshot of the velocity field for the initial and an intermediate time step are shown. In it, the elastic behaviour of the fluid is pointed out. In this case, some minor loss of symmetry is observed. This reveals that the number of particles associated to each node is not large enough. We verified that by increasing the number of molecules the symmetry is enhanced.

There is another reason for this asymmetry of a more geometrical nature. We employed three-points integration on the triangles. As is well-known, Delaunay triangulation for four nodes lying on the same circle is not unique, and the use of a very regular node distribution at the early stages of the simulation can lead to non-symmetric triangulations at both sides of the symmetry axes. This can be minimised by using Stabilised Conforming Nodal Integration techniques [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF].

Die swelling for the Doi-Edwards model

The same geometry used in Section 6.2 is again employed for the simulation of die swelling of an entangled polymer modelled with the Doi-Edwards theory presented above. In this case, a time step of 0.007s was employed with λ H = 1.0, and only 100 molecules were attached to each node to alleviate the computational cost involved in this kind of micro-macro simulations.

Four snapshots of the velocity field are shown in Fig. 7. Again, elastic effects are notorious after the outlet of the channel. A minor loss of symmetry in the flow is noticed due mainly to the statistical noise. Despite the very low number of realisations per node, the statistical noise remains surprisingly low.

Conclusions

A method for the simulation of complex models arising from kinetic theory of fluids has been presented. This method is based on the combination of two different ingredients. On one hand, the use of a CONNFFESSIT-like strategy, in which benefit is obtained from the equivalence of the Fokker-Planck equation and an Itô's stochastic differential equation. Secondly, an updated Lagrangian approach by employing the Natural Element Method. This meshless strategy allows to associate the particles carrying the realizations of the stochastic process to the nodes. They will consequently move with material velocity, avoiding the need of complex particle tracking algorithms and the related numerical problems (lack of particles in certain regions of the domain, etc.) The technique is used in conjunction with shape constructors to track the evolution of the free surface.

A comparison with existing results for FENE models has been presented. It shows good accuracy, with respect to results existing in the literature for homogeneous flows. Some examples of free-surface die swelling flows with both FENE and Doi-Edwards models have also been presented. Although no closed-form results exist for these problems, the results obtained show a priori a promising potential of the proposed technique in the prediction of complex, non-homogeneous flows in the presence of free surfaces.

The main drawback of the technique is, maybe, the computational cost. In order to keep the number of degrees of freedom reasonably low, a moderate number of realisations of the stochastic process should be accomplished. The price to pay is an obvious increase in the numercial noise of the results, that can eventually lead to some minor loss of symmetry in certain simulations. The FENE model has demonstrated to be more sensitive to this effect than the Doi-Edwards model. We believe, however, that the presented technique constitute a valuable alternative to traditional finite elements, since it avoids both the mesh distortion effects associated to Eulerian and ALE approaches, and the numerical diffusion associated to remeshing. In any case the stochastic integration that constitutes the most expensive computational operation can be carried out in parallel using appropriate computing platforms.

Figure 1 :

 1 Figure 1: Delaunay triangulation and Voronoi diagram of a cloud of points.

Figure 2 :

 2 Figure 2: Definition of the Natural Neighbour coordinates of a point x.

Figure 3 :

 3 Figure 3: Evolution of the family of α-shapes of a cloud of points representing a wave breaking on a beach. Shapes S 0 (a), S 0.5 (b), S 1.0 (c), S 2.0 (d), S 3.0 (e) and S ∞ (f) are depicted. 8

Step 1 :• Step 2 :

 12 Solve the macro equations via FEM and obtain the velocity field v Using v, compute the path of each molecule • Step 3: Integrate the S.D.E. governing the micro level • Step 4: Update the viscoelastic stress by means of averaging

Step 3 :

 3 For each node, integrate the S.D.E. and update configuration • Step 4: Update the viscoelastic stress by means of averaging

Figure 4 :

 4 Figure 4: Evolution of non-dimensional viscosity (a), first stress normal coefficient (b) and second normal stress coefficient (c) in a FENE fluid undergoing the startup of shear flow.

Figure 5 :

 5 Figure 5: Geometry and initial configuration of nodes in the channel.

Figure 6 :

 6 Figure 6: Velocity field in the FENE die swelling simulation for two different time steps.

Figure 7 :

 7 Figure 7: Evolution of the velocity field in the die swelling flow of reptating polymer.

Even in three-dimensional spaces, it is common to refer to the Delaunay tetrahedralisation with the word triangulation in the vast majority of the literature

The α-shapes-based Natural Element MethodThe identification of the free surface in an updated Lagrangian flow simulation deserves some comments. In many prior works, location of boundary nodes is performed by flagging coincident element faces[START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF], for instance. Once the updating of nodal positions has been performed, a recursive check must be done in order to find overlapping boundary segments. In three di-

Acknowledgements

This work has been partially supported by the Spanish Ministry of Science and Innovation, through grant number CICYT-DPI2008-00918. This support is gratefully acknowledged.