Saumen Dey

Khalid Belhajjame

David Koop

Tianhong Song

Paolo Missier

Bertram Ludäscher

UP & DOWN: Improving Provenance Precision by Combining Workflowand Trace-Level Information

Workflow-level provenance declarations can improve the precision of coarse provenance traces by reducing the number of "false" dependencies (not every output of a step depends on every input). Conversely, fine-grained execution provenance can be used to improve the precision of input-output dependencies of workflow actors. We present a new logic-based approach for improving provenance precision by combining downward and upward inference, i.e., from workflows to traces and vice versa.

Introduction

Many scientific workflow systems have been instrumented to capture workflow execution events as provenance. Typically, such events record the computational steps that were invoked as part of a workflow execution as well as the data that were input to and output by each step. These recorded events can be used to trace data provenance by identifying the input values that contributed to an output data value generated as a result of the workflow execution. In practice, however, this application is limited by the black box nature of the actorsthe computational modules that implement the steps of the workflow.

To illustrate this, consider the actor A shown in Fig. 1. Assume that an execution of A uses the input values v I 1 and v I 2 for the input ports I 1 and I 2 , respectively, and output ports O 1 and O 2 produce values v O 1 and v O 2 , respectively. Because the actor is a black box we cannot assert whether the value v O 1 was derived from (i) v I 1 , (ii) v I 2 , (iii) both v I 1 and v I 2 , or (iv) none of the inputs. Similarly, we can not assert how the value v O 2 was derived from v I 1 and v I 2 . By derivation, we refer to the transformation of a data value into another, an update of a data value resulting in a new one, or the construction of a new data value based on a pre-existing one [START_REF] Moreau | PROV-DM: The PROV data model[END_REF]. All we can safely conclude is that the invocation of the actor used the input values v I 1 and v I 2 and generated the values v O 1 and v O 2 . While useful, this information is not sufficient to trace fine grained dependencies of data values produced by workflow executions. For the actor A shown * University of California, {scdey,thsong,ludaesch}@ucdavis.edu † Université Paris-Dauphine, Khalid.Belhajjame@dauphine.fr ‡ New York University, dakoop@nyu.edu § Newcastle University, Paolo.Missier@ncl.ac.uk

I 1 I 2 O 1 O 2
(a)

I 1 I 2 O 1 O 2 (b) I 1 I 2 O 1 O 2 (c)
Figure 2: Three examples of internal port dependencies of the actor A shown in Fig. 1. There could be 13 more such dependencies.

in Fig. 1, there could be 16 possible ways output ports are dependent on the input ports, three of which are shown in Fig. 2. If we determine that one of the input values was incorrect, e.g., due to a malfunctioning sensor, we would potentially need to invalidate many useful results than if we knew that the input was only used for an insignificant output. This is one of many reasons showing the importance of finding fine-grained dependencies.

In general, if there are n input ports and m output ports for an actor, then there are 2 mn possible internal port dependencies. With k such actors in a workflow, there are 2 kmn possible internal port dependency models of the workflow, i.e., possible models of the workflow.

In this paper, we propose a framework that takes a workflow specification and a set of provenance information to generate all possible data dependency models. By dependency models, we mean a graph in which the nodes represent the input and output ports of the actors, and the edges specify derivation dependencies between the ports. Notice that the framework generates multiple dependency models, all of which are possible. Fig. 2 shows three possible dependency models of the workflow with a single actor A shown in Fig. 1. However, given the execution trace of a workflow execution, only one of the possible data dependency models reflects the true dependencies between the ports, whereas the remaining data dependency models are false positives. Figure 3 depicts the overall framework of our proposal towards reducing the number of false positive data dependency models. During the design of a workflow, the designer knows some internal details of the actors, e.g., (i) data on an output port can only be produced only after using data from a specific subset of input ports, (ii) two outputs are generated using exact set of inputs, etc. In step 1 , we assume that our workflow model captures these additional design level information and analyze them to generate additional dependencies. In 2 , we analyze the provenance graph structure and infer further dependencies by computing the transitive closure of data dependencies. Given a set of provenances either by (i) collecting actual workflow execution traces, or (ii) probing the input values, executing the workflow, and collecting provenances from the execution traces, in 3 we analyze input and output data values towards understanding the inputoutput dependencies. Steps 1 , 2 , and 3 generate additional dependency information, which this framework encodes as constraints and uses them while generating all the possible models of the workflow in step 4 . Each of these models is an "improved" version of the given workflow, i.e., they have more dependency information, and these models are in turn used in step 5 to validate a provenance graph of an unknown origin (i.e., workflow) and improve it by injecting more dependencies. This improved provenance then can be used to have better result analysis, debugging, etc.

Running Example. We use the Climate Data Collector workflow to showcase the features of the proposed framework. We assume a set of sensors that collect temperature, pressure, and other climate observations including various flags. Given a sensor id, the ReadSensor actor reads from the sensors and returns temperature, pressure, and three flags, flagA, flagB, and flagC. The second actor, SensorLogic, takes these three flags and returns two flags weatherCode and temperatureCode. The flag weatherCode specifies if the weather was normal or not and temperatureCode indicates if the temperature was read in Celsius or in Fahrenheit. Another actor, ConvertToKelvin, uses the temperatureCode flag and the temperature reading to convert the temperature to Kelvin. A final actor RangeCalculation accepts the pressure and temperature and computes a range. This workflow, defined in VisTrails, is shown in Figure 4.

If we consider all these actors as black boxes, ReadSensor, SensorLogic, ConvertToKelvin, and RangeCalculation actors individually have 32, 64, 4, and 4 input-output dependency models and thus, the workflow as a whole has 32768 possible dependency models. These many models makes it impossible to make any specific claims about dependencies, but if we use provenance information, we can reduce the number of possibilities.

For example, all the outputs of ReadSensor depend on a single input and there is no output if the sensor is not available. These information reduces 32 possible dependency models to exactly one which makes result analysis easier. Unfortunately, neither the workflow specification nor the provenance clearly capture these dependencies. We use this example workflow and show how we can capture additional dependency information, how we can use this information, and how we can improve the possible dependency models in Section 4.

The paper is organized as follows. We start by stating the assumptions needed for workflow model and provenance model used in our approach in Section 2. We then explain our proposed framework and the logical architecture in Section 3. In Section 4, we discuss our implementation using the example workflow and then describe related work in Section 5. We conclude the paper and discuss ongoing work in Section 6.

Models

Workflow Model. We base our workflow model on [START_REF] Cuevas-Vicenttín | Modeling and querying scientific workflow provenance in the dopm[END_REF][START_REF] Cuevas-Vicenttín | The PBase scientific workflow provenance repository[END_REF]. A workflow W = (V w , E w ) is a directed graph whose nodes V w = A ∪ P ∪C are actors A, ports P, and channels

Proposed Framework and Architecture

This framework takes a workflow specification W and a set of provenance graphs Gs and infers the internal (i.e., internal to an actor) port dependencies by analyzing W and Gs. It then encodes all these inferred dependencies as constraints, which are used to reduce the number of possible workflow models. It uses the Generate-and-Test pattern from Answer Set Programming (ASP) [START_REF] Lifschitz | What is answer set programming[END_REF], i.e., we use sets of stable models [START_REF] Gelfond | The stable model semantics for logic programming[END_REF] to represent all possible models that are consistent with the generated constraints.

Let us consider the snippet of ASP program we have used in this framework.

pd(X,Y) v pnd(X,Y):-in(X,A), out(A,Y).

Here, pd(X,Y ) relation means that Y depends on X and npd(X,Y ) relation means that Y does not depend on X. If there are n in(X, A) relations and m out(A,Y ) relations in W , then there would be 2 mn dependency models for the actor A. In one extreme possible model, we would have all the pd(X,Y ) as "true" and in another extreme possible model, we would have all the npd(X,Y ) as "true". In all possible models, we would have some of pd(X,Y ) as "true" and some npd(X,Y ) as "true". Now, this framework applies the constraints in the following way: if a possible model has one npd(X,Y ) as "true" and for the same (X,Y ) pair a constraint has been observed to be "true" as discussed in Section 1 then its a contradiction as the constraint specifies that node Y depends on node X. The framework applies the negation as failure principle and excludes this model.

The proposed framework, as shown in Fig. 5, has several components including the Constraint Generator, Model Generator, Model Reducer, and Trace Validator. The Constraint Generator takes a workflow W and a set of provenance traces Gs, and performs three analysis tasks to generate constraints using the techniques we present next.

We analyze provenance graphs and capture such inferred data dependencies and use them to reduce the number of models. More specifically, in our provenance model, we capture used(I, D), genBy(D, I), and derBy(D 1 , D 2 ) edges. The used(I 2 , D), and genBy(D, I 1 ) edges along with the mappings invoc(I 1 , A 1 ), and invoc(I 2 , A 2 ) unambiguously states that the actor A 2 depends on actor A 1 . However, edges used(I, D 1 ), and genBy(D 2 , I) are unable to unambiguously state that data artifact D 2 depends on data artifact D 1 . The derBy(D 1 , D 2 ) edges describe dependencies among data artifacts and thus these edges are able to infer internal port dependencies of an actor. This framework analyzes these derBy/2 (i.e., derBy is a two arity relation) relations and infers further data dependencies.

We use the provenance graph shown in Fig. 6 to describe how this framework analyzes these derBy/2) relations. This provenance graph has seven data artifacts T though Z and two derBy/2 edges, derBy(Z,Y ) and derBy(Z, T ) . Considering the derBy(Z, T ) edge we infer that Z depends on W . Thus, given the derBy(Z,Y ) and derBy(Z, T ) edges, we conclude that output port associated with Z depends on input ports associated with W and Y , which reduces the number of possible dependency models. More generally, given provenance traces Gs, the derBy/2 relations in all the Gs are analyzed. A derBy edge could specify the dependency of (i) an output and an input of some actor or (ii) an output of one actor and an input of another actor. The framework infers these constraints using an algorithm, which is represented by the datalog rules 4 and 5 as shown below. Based on these two inferences, we get that port P2 depends on P1, which is a workflow specification level information by analyzing the derBy/2 relation from the provenance graph. Note that (i) provide exact dependency, but (ii) provides an over estimate. Input port dependencies can also be inferred by analyzing the input and output data values * . The framework computes the constraints using an algorithm, which is represented using rule 6 shown below. The Model Generator computes all possible dependency models in which output ports may depend on input ports for an actor and removes the models that contradict the constraints generated by the Constraint Generator as discussed in the beginning of this section. This component returns only those models that have no contradictions. The Model Reducer takes all the possible workflow models generated by Model Generator and combines them into one workflow model. This workflow model is an improved specification of the input workflow * In this work, we assumed a simple model for input and output data value analysis. In this model, an output port of an actor depends on an input port if the values on the output port changes while changing values of an input and keeping the values other inputs same. In our future work, we plan to investigate further complexities of an actor. model W , in which some of overestimates of the internal port dependencies have been removed. In Provenance Validator, the framework uses this improved workflow specification (i) to validate the provenance graph of unknown origin, (ii) to improve provenance graph by removing the inconsistencies.

T U A 2 A 1 A 4 V W A 3 X Z Y

Prototypical Implementation

In this section, we describe a prototypical implementation of the proposed framework. Fig. 7 shows the Climate Data Collector workflow represented using this framework; actors are represented as rectangles and have ports "p01" through "p18". The edges represent channels which connect output ports of one actor to input ports of another actor.

Fig. 8 shows that there are 32 possible dependency models for the SensorLogic actor. Without additional information, the framework will generate all 32 models. Similarly, other actors have many possible dependency models (as discussed in Section 1), and are not displayed because they mirror the SensorLogic actor.

The framework analyzes the workflow specification, provenance traces, and user-specified constraints as discussed in Section 1, and generates the constraints as shown in Fig. 9 which the framework applies while generating the possible dependency models. Fig. 9(a) shows the constraints generated based on the additional specifications † provided by the workflow designer. Fig. 9(b) shows all the constraints generated by the framework by analyzing the input-output value dependencies. Fig. 9(c) and Fig. 9(d) show all the constraints generated by the framework by analyzing the derBy/2 edges from the provenance graphs. Fig. 10 shows the improved workflow specification based on the given workflow, provenance Using this example, we have demonstrated that the proposed framework can improve the workflow specification, and this information can be used to later validate provenance graphs with unknown origin.

Related Work

The problem of defining or inferring data dependencies in scientific workflows has been investigated by a handful of researchers. For example, Bowers et al. proposed a declarative language for specifying fine-grained dependencies at the level of the workflow definition, which are propagated and applied to workflow trace events produced as a result of the workflow execution [START_REF] Bowers | Declarative rules for inferring fine-grained data provenance from scientific workflow execution traces[END_REF]. Their approach is complementary to ours, and we adopt a similar language to encode the dependencies that are dervied from workflow specifications and their corresponding trace events. However, that approach did not tackle the problem of inferring data dependencies. Ghoshal et al. investigated the use of static analysis techniques to derive dependencies (or what they term mappings) between the inputs and outputs of an actor [START_REF] Ghoshal | Static compiler analysis for workflow provenance[END_REF]. However, while they assume that the source code of the program implementing an actor is available, we tackle the problem of deriving data dependencies when the source code for the actors is not avaiable.

Garijo et al. proposed a framework where they identified data transformation and manipulation patterns (which they term motifs) that are commonly found in scientific workflows [START_REF] Garijo | Common motifs in scientific workflows: An empirical analysis[END_REF]. They distinguish two types of motifs: data intensive activities that are observed in workflows (data-oriented motifs), and different manners in which activities are implemented within workflows (workflow-oriented motifs). The first type is relevant to our work. However, in our framework, we are particularly interested in identifying if there is a dependency between an input and an output of an actor rather than the kind of data manipulation performed by the actor.

The above work and others, e.g., [START_REF] Alper | Enhancing and abstracting scientific workflow provenance for data publishing[END_REF], [START_REF] Sethi | Re-using workflow fragments across multiple data domains[END_REF], are related because they aim to understanding the kind of manipulation carried out by workflows. However, none tackle the problem of identifying fine-grained dependencies between input and output port for black box actors. In our prior work [START_REF] Dey | Provenance analyzer: exploring provenance semantics with logic rules[END_REF], we used the same Generate-and-test ASP based approach towards finding the possible orders of events. However, in this work we focus on finding the fine-grained dependencies using the same Generate-andtest ASP based approach.

Discussion

In this paper, we have presented a framework for inferring fine grained dependencies between the inputs and output ports of actors. We have described a probing based methods for inferring such dependencies. In doing so, we made the assumption that dependencies are static, in the sense that they hold across all invocations of the actors. In practice, however, dependencies may be dynamic, in the sense they may change across invocations of the actor. For example, the VisTrails system [START_REF] Freire | The Architecture of Open Source Applications, chapter VisTrails[END_REF] provides a conditional actor (If) that uses modified dataflow logic to execute only one of two upstream workflows to generate an output based on a boolean input. There are also instances where two ports are provided to make the specification of an input possible in different formats, meaning only one port's value will eventually be used but is dependent on which are provided.

The framework that we presented in this paper needs to be refined in order to cater for the identification of dynamic dependencies. In particular, the user should be able to understand the cases (conditions under which) a given output is likely to depend or not on a given input. Note also, that so far, we treated input-output port equally. A classification of dependencies is needed to provide the user (workflow designer or provenance user) with better understading on the kind of relationship between the input and output ports. For example, distinguishing control-flow dependencies from data-flow dependencies. The later can be further classified to specify the kind of contribution a given input value had in the construction of the output value.

2 Figure 1 :

 21 Figure 1: An actor A with two input ports I 1 , and I 2 and two output ports O 1 and O 2 . Do the outputs depend on none, one, or both of the inputs?

Figure 3 :

 3 Figure 3: Components of our solution for inferring dependencies.

Figure 4 :

 4 Figure 4: Climate Data Collector workflow, an example workflow that processes climate data read from a sensor.

Figure 5 :

 5 Figure 5: Logical architecture of how components interact to improve a workflow specification.

( 1 )

 1 ddep(X,Y) :-used(I,X),genBy(Y,I). (2) ddep * (X,Y) :-ddep(X,Y). (3) ddep * (X,Y) :-ddep(X,Z),ddep * (Z,Y).

Figure 6 :

 6 Figure 6: A provenance graph with two derived-by edges derBy(Z,Y) and derBy(Z,T). Circles and rectangles represent data artifacts and processes, respectively.

( 4 )

 4 pdep(P1,P2) :-derBy(Y1,X1),data(X1,P1,A,R,_), data(Y1,P2,A,R,_). (5) pdep(P1,P2) :-derBy(Y1,X1),data(X1,P1,A,R,_), data(D,P2,A,R,_),ddep * (D,Y1), ddep * (X1,D).

( 6 )

 6 pdep(IP,OP) :data(D1,IP,A,R1,IV1),data(D1,IP,A,R2,IV2), data(D2,OP,A,R1,OV1),data(D2,OP,A,R2,OV2), ¬ R1=R2, ¬ IV1=IV2, ¬ OV1=OV2.

Figure 7 :

 7 Figure 7: Climate Data Collector workflow represented by this framework. There are 32768 possible dependency models for this workflow.

Figure 8 :

 8 Figure 8: Possible dependency models for the Sensor-Logic actor. There are 32 possible dependency models for this actor alone.

Figure 9 :

 9 Figure 9: All of the constraints that are generated by the framework.

Figure 10 :

 10 Figure 10: Improved workflow specification of the Climate Data Collector workflow after including one possible dependency model. There are be 6 such possible models, and this framework generates all of them.

† The definition of this specification is available in[START_REF] Dey | Datalog as a lingua franca for provenance querying and reasoning[END_REF].

Acknowledgments. Supported in part by NSF ACI-0830944 and IIS-1118088.


Angelo Simone
email: a.simone@tudelft.nl

Partition of unity-based discontinuous finite elements: GFEM, PUFEM, XFEM

Keywords: discontinuités en déplacement, GFEM, PUFEM, XFEM displacement discontinuities, GFEM, PUFEM, XFEM

Dans cet article, nous présentons les principes élémentaires des éléments finis discontinus basés sur la méthode de la partition de l'unité en montrant leur relation avec la méthode des éléments finis généralisée. Les aspects liés à l'implantation dans un code ainsi que la simplicité de cette méthode sont illustrés sur un exemple unidimensionnel. La capacité de cette approche à représenter des discontinuités en déplacement indépendamment du maillage élément fini est montrée sur la propagation d'une fissure dans un milieu élastique. Nous montrons aussi certaines limitations de cette méthode quand elle est utilisée avec l'approche "dummy stiffness". ABSTRACT. In this paper we review some basic notions of partition of unity-based discontinuous finite elements showing their relation to the Generalized Finite Element Method. A minimal one-dimensional example illustrates some of the issues related to the computer implementation of the method and highlights the relative simplicity of the approach. The ability of the approach in describing displacement discontinuities independently of the finite element mesh is shown in a classical crack propagation problem in an elastic medium. We also illustrate some limitations of this method when used in conjunction with the dummy stiffness approach.

Introduction

Displacement discontinuities are incorporated in standard finite elements for the numerical treatment of geometrical discontinuities, like e.g. rock joints [START_REF] Goodman | A model for the mechanics of jointed rocks[END_REF], or for the description of problems with evolving boundaries, like e.g. cracks in quasi-brittle materials [START_REF] Saouma | A review of fracture mechanics applied to concrete dams[END_REF] or delamination in composite materials (Alfano et al., 2001). The incorporation of a discontinuity in the displacement field is conventionally pursued by using interface elements. With interface elements, a discontinuous displacement field is described by the relative displacement of a double set of nodes at inter-element boundaries. This requires special mesh generators and, in the case of evolving boundaries, remeshing procedures [START_REF] Bittencourt | Quasi-automatic simulation of crack propagation for 2D LEFM problems[END_REF]. Alternatively, a discontinuous displacement field can be described by means of embedded discontinuity elements in which the displacement discontinuity, described by using the Heaviside function, is incorporated into the finite-element formulation as an incompatible strain mode and can pass through solids elements [START_REF] Simo | An analysis of strong discontinuities induced by strainsoftening in rate-independent inelastic solids[END_REF][START_REF] Jirásek | Comparative study on finite elements with embedded cracks[END_REF].

A different approach employs the partition of unity property of finite-element shape functions (the sum of the shape functions must equal unity at each spatial point). Within this approach, the standard approximation basis is enriched locally with special functions. This enrichment results in extra degrees of freedom for the nodes in the domain subjected to the enrichment, without modification of the mesh topology (Babuška et al., 1994;[START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF][START_REF] Babuška | The partition of unity finite element method[END_REF]Duarte et al., 1996b;Duarte et al., 1996a;[START_REF] Oden | A new cloud-based hp finite element method[END_REF][START_REF] Oden | Clouds, cracks and FEM's[END_REF]. The effect of a displacement discontinuity is described by enriching the standard finite-element polynomial basis with a step function [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. The outcome of this approach is a class of elements, known in the literature under the names "generalized finite element method" (GFEM) [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF], "partition of unity finite element method" (PUFEM) [START_REF] Gasser | Modeling 3D crack propagation in unreinforced concrete using PUFEM[END_REF] and "extended finite element method" (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], which is kinematically equivalent to the class of conventional interface elements, the key difference being the possibility of arbitrarily locating the discontinuity within the domain of an element. When PU-based discontinuous elements are used, the interface response is described by an additional set of global degrees of freedom and by a constitutive law at the discontinuity.

Beyond the classical finite element method: Exploiting partitions of unity

In this section, we recall some basic notions on the partition of unity paradigm and its relation to finite element methods. A partition of unity is defined as a collection of global functions N I (x) whose value sums up to unity at each point x in the solution domain V of a boundary value problem:

n I=1 N I (x) = 1 ∀x ∈ V . [1]
where n is the number of nodes of the discretization. In the finite element method (FEM) [1] holds by construction for any Lagrangian FEM global basis (i.e. shape function). Using linear partitions of unity N I as basis functions in a FEM framework [START_REF] Taylor | A hierarchical finite element method based on the partition of unity[END_REF], a scalar field q(x) can be approximated by

q h (x) = n I=1 N I (x)q I regular interpolation + n I=1 N I (x) m j=1 E j (x)q Ij enrichment , [2]
where the subscript "h" denotes the approximate field, qI (termed "regular" or "standard" degrees of freedom (dofs)) are nodal values related to the basis N I , and qIj (termed "enrichment" or "extra" degrees of freedom) are the global nodal parameters connected to the basis E j pertinent to the enrichment (m is the number of terms in E j ).

The terms "regular" and "enrichment" make reference to the fact that the "regular" interpolation field is considered as the background field upon which the "enrichment" is added. To avoid linear dependency, the order of any polynomial terms in the enrichment basis must be greater than the order of the partitions of unity N I . This is, however, not a sufficient condition; algorithms to deal with the linear dependent case are described in [START_REF] Duarte | Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems[END_REF].

Shape functions associated with [2] are termed GFEM shape functions. These shape functions are built from the product of a partition of unity N I and enrichment functions E j :

N I × {1, E j }.
[3]

In the GFEM, the partition of unity is in general provided by linear Lagrangian FEM shape functions. The partition of unity property of the functions N I implies that

n I=1 N I (x)E j (x) = E j (x). [4]
Consequently, any enrichment function E j can be represented exactly through linear combinations of GFEM shape functions, and, if the enrichment functions E j can approximate the solution of a boundary value problem, the corresponding GFEM shape functions also will.

In finite element notation, the approximation for a vector field q(x) of order l over an n-node enriched element, with all the nodes enriched, can be written as

q h (x) = N (x) q + N (x)N E (x) q, [5]
where N is a l × (l × n) matrix containing the standard finite element shape functions, N E is a (l×n)×(l×m×n) matrix containing the extra basis terms for the enrichment, q is a (l × n) × 1 vector containing standard degrees of freedom and q is a (l × m × n) × 1 vector containing the extra degrees of freedom. The number of extra degrees of freedom per node is equal to the number of terms in the basis for the enrichment multiplied by the number of nodal unknowns. In standard finite elements, the matrix N E is empty.

As noted in [START_REF] Oden | A new cloud-based hp finite element method[END_REF] and [START_REF] Taylor | A hierarchical finite element method based on the partition of unity[END_REF], the approach of Equation [5] allows the enrichment to be performed from node to node in a mesh by activating the extra degrees of freedom q when needed (a hierarchical finite element formulation based on the partition of unity method).

Enrichment functions

The essence of the generalized FEM is condensed in [2]. This special approximation allows to paste a-priori knowledge about the behavior of the solution of a problem directly into the shape functions. Using this feature, polynomial functions, discontinuous functions and even singular functions have been exploited as enrichments. These customized enrichment functions have been used to model local features like cracks, voids, inclusions, microstructures etc.

Figure 1 shows the construction of a generalized FEM shape function as the product of a partition of unity and an enrichment function for two different enrichment functions. Next, we concentrate on the discontinuous enrichment.

Kinematics for a body crossed by a discontinuity

To incorporate a discontinuity in the displacement field it is necessary to characterize the kinematics of a body crossed by a discontinuity. A body V bounded by S and crossed by a single discontinuity surface S d is considered (see Figure 2). The internal discontinuity surface S d divides the body into two sub-domains, V + and V - (V = V + ∪ V -). Displacements ū are prescribed on S u , while tractions T are prescribed on S t . The displacement field can be decomposed as [6] where H(x) is the Heaviside function centered at the discontinuity surface

u(x, t) = û(x, t) + H(x)ũ(x, t) if x / ∈ S d ,
S d (H = 1 if x ∈ V + , H = 0 if x ∈ V -
) and û and ũ are continuous functions on V . The discontinuity is introduced by the Heaviside function H at the discontinuity surface S d and the value of the displacement jump u across the discontinuity surface is given by ũ at S d . In the geometrically linear case, the strain field in V is computed as the symmetric part of the gradient of the displacement field and reads [7] where (•) s refers to the symmetric part of (•). The strain field is defined everywhere in the body except at the discontinuity surface.

ε (x, t) = ∇ s û (x, t) + H (x) ∇ s ũ (x, t) if x / ∈ S d ,
S u S u V - n t S d n d V + S t Figure 2. Body V crossed by a discontinuity S d

Describing a displacement discontinuity

In describing a displacement discontinuity, and considering [5] now rewritten for convenience as

u h = N a + N N E b, [8]
if the standard displacement field is interpolated by the regular interpolation N a, the displacement jump can be described through the enrichment N N E b = HN Hb [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], where H is the scalar valued Heaviside function (m = 1) and H is a diagonal matrix identifying through a 1/0 switch which degrees of freedom to enrich (H is the identity matrix if all the degrees of freedom of the element are enriched). This results in the (l × n) × 1 vector of extra degrees of freedom b (same dimensions as the vector of standard degrees of freedom a) and in the (l × n) × (l × n) matrix HH. The enrichment concerns only nodes whose support is crossed by a discontinuity. For nodes whose support is not crossed by a discontinuity, the basis responsible for the enrichment is empty since the Heaviside function is a constant function over their supports and can be neglected. Since element stiffness matrices for elements with extra degrees of freedom are assembled only for active extra degrees of freedom, the standard matrix N can be used in place of N H. Hence, the discretized format of [6] reads, for points located in an element with enriched nodes,

u h = N a + HN b, [9]
where the global nodal degrees of freedom a and b represent, in the arrangement of [9], the total displacement field. The displacement jump across the discontinuity S d is given by

u h = N b ⏐ ⏐ ⏐ S d
.

[10]

For points located in an element without enriched nodes, the Heaviside function is a constant function over their supports and therefore it is not considered. Consequently, since there is no enhancement, the standard finite-element interpolation u h = N a is retrieved.

It is worth noting that in this case the GFEM shape functions [3] are built from the product of a partition of unity N I and the step function H as shown in Figure 1b:

N I × {1, H}.
[11]

The case of multiple discontinuities is easily dealt with by adding the corresponding enrichment functions (a detailed description of the procedure in case of intersecting discontinuities is reported in [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF] and [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF]).

Near-tip enrichment

So far, we have considered the enrichment across a crack surface. However, with reference to the crack geometry depicted in Figure 3 care should be taken to properly describe the singularities at the crack tip in case of a linear elastic sample. As proposed by [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF], a set of near-tip functions that spans the exact asymptotic displacement field near a crack-tip in a two dimensional linear elastic body subjected to general mixed-mode loading can be used as enrichment around the crack tip. One such a set is given by the functions

E (r, θ) = √ r sin θ 2 , √ r cos θ 2 , √ r sin θ 2 sin θ, √ r cos θ 2 sin θ , [12]
where (r, θ) are the local polar coordinates at the crack as shown in Figure 3.

With these functions as enrichment at the crack tip, the complete displacement field is expressed as [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF])

u h = i∈I a i N i continuous + j∈J b j N j H discontinuous + k∈K N k 4 l=1 c kl E l tip , [ 13
]
where I is the set of standard nodes, J is the set of circled nodes and K is the set of squared nodes reported in Figure 3. Note that the use of the functions reported in [12] is limited to two-dimensional linear elastic fracture mechanics and should not be used in other situations. Also, the inclusion of the near-tip field in [12] is an important advantage over other approaches in which only the discontinuity in the displacement field is considered (cf. embedded discontinuities models [START_REF] Jirásek | Comparative study on finite elements with embedded cracks[END_REF]).

Hp discretization

Near-tip enrichment functions are not always available. Indeed, the stress state in the neighborhood of a crack front is not well known in three-dimensions, and analyti-cal expansions are available only for particular crack geometries [START_REF] Broberg | Cracks and Fracture[END_REF][START_REF] Mittelstedt | Semi-analytical computation of 3D stress singularities in linear elasticity[END_REF]. Complex crack front geometries, curved crack surfaces or the intersection of the crack surface with the boundary, create complex stress distributions that are, in general, not amenable to closed form expansions. Similarly, non-linear fracture models generate stress states that are not amenable to analytical expansions in the neighborhood of a process zone. Two additional difficulties with near crack tip enrichment functions are their numerical integration and the need to enrich several layers of elements around the crack front in order to achieve optimal convergence rates [START_REF] Laborde | High-order extended finite element method for cracked domains[END_REF][START_REF] Béchet | Improved implementation and robustness study of the X-FEM for stress analysis around cracks[END_REF][START_REF] Xiao | Improving the accuracy of XFEM crack tip fields using higher order quadrature and statically admissible stress recovery[END_REF]. Integration errors may lead to wrong conclusions in a convergence analysis since they may cancel discretization errors when computing, for example, the strain energy.

A strategy to cope with the lack of general near-tip enrichments, consists in resolving the non-smooth displacement field near a crack tip/front by local mesh refinements combined with high-order approximations [START_REF] Szabó | Mesh design for the p-version of the finite element method[END_REF]) (hp discretization). This approach avoids the difficulties and limitations related to near crack front enrichment functions while being able to deliver high convergence rates as demonstrated in Duarte et al. (2007). The construction of hp discretizations in the neighborhood of a crack front is not difficult in the generalized FEM since crack surfaces are allowed to cut elements. The mesh refinement can, therefore, be performed as if there were no cracks in the domain.

High-order approximations

At variance with the established enrichment of the displacement field with discontinuous functions and the use of linear approximations proposed by [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] in the so-called XFEM (see [13]), approximation orders higher than one have been exploited in various forms by [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF], [START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF], Chessa et al. (2003), [START_REF] Mariani | Extended finite element method for quasi-brittle fracture[END_REF], and [START_REF] Laborde | High-order extended finite element method for cracked domains[END_REF].

The performance of the method can be improved considerably by consistently using high-order hierarchical GFEM approximations based on a linear finite element partition of unity. In this situation, all nodal degrees of freedom are defined at element vertex nodes which facilitates mesh generation, especially in three-dimensional computations. This is in contrast with the use of high-order non-hierarchical Lagrangian partition of unity functions with piece-wise constant enrichment functions [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF][START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF][START_REF] Chessa | On the construction of blending elements for local partition of unity enriched finite elements[END_REF][START_REF] Laborde | High-order extended finite element method for cracked domains[END_REF]. This approach leads to nodes on element edges, faces and interior. Hierarchical approximations are also ideal for the construction of non-uniform polynomial approximations which can be very effective for some classes of problems [START_REF] Duarte | Mesh-independent directional p-enrichment using the generalized finite element method[END_REF]. Other advantages of hierarchical GFEM, such as the non-zero structure of the resulting stiffness and mass matrices, are discussed in [START_REF] Taylor | A hierarchical finite element method based on the partition of unity[END_REF], [START_REF] Oden | A new cloud-based hp finite element method[END_REF], and [START_REF] Duarte | Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems[END_REF]. The issue of linear dependence of hierarchical GFEM enrichments is also discussed in detail in [START_REF] Duarte | Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems[END_REF].

Below, we list a few examples of GFEM approximations for single cracks in twodimensions. Hereafter, N I I = 1, ..., n are standard Lagrangian shape functions, x I = (x I , y I ) are the coordinates of the node I, h I is the diameter of the largest finite element sharing the node I, n is the number of nodes in the mesh, p and q are used to denote the polynomial degree of ûh and ũh , respectively. 1) p = 1, q = 1: the shape functions at a node x I are given by

N I × {1, H} [14]
which are the same generalized FE shape functions in [11] i.e. we recover the simplest possible case, that of a piecewise linear approximation of u. These shape function will be used in the remainder of the paper.

2) p = 2, q = 0: this case corresponds to the generalized FE shape functions introduced in [START_REF] Oden | A new cloud-based hp finite element method[END_REF] and [START_REF] Duarte | Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems[END_REF]. The shape functions at a node x I are given by

N I × 1, (x -x I ) h I , (y -y I ) h I .
[15]

These shape functions, of course, can not represent the discontinuous components of the solution ũh .

3) p = 2, q = 2: the shape functions at a node x I are given by

N I × 1, (x -x I ) h I , (y -y I ) h I , H, H (x -x I ) h I , H (y -y I ) h I .
[16] 4) p = 2, q = 1: the shape functions at a node x I are given by

N I × 1, (x -x I ) h I , (y -y I ) h I , H . [17]
These shape functions are equivalent to those proposed by [START_REF] Stazi | An extended finite element method with higher-order elements for curved cracks[END_REF].

A study of the performance of the above shape functions can be found in Duarte et al. (2007).

Governing equations

Having characterized the displacement field for a body crossed by a discontinuity, we now briefly recall the governing equations for such a body. The equilibrium equations and boundary conditions for the body V without body forces depicted in Figure 2 can be summarized by

∇•σ = 0 in V, [18a] σn = t on S t , [18b] σn d = t on S d , [18c]
where σ is the Cauchy stress tensor, n is the outward unit normal to the body and n d is the inward unit normal to V + on S d . Equation [18c] represents traction continuity at the discontinuity surface S d . The strong form is completed by the boundary conditions [19b] where ū is a prescribed displacement. Equation [19a] is a standard essential boundary condition while [19b] has been imposed to simplify the finite-element implementation [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF]. The constitutive relationships for bulk and discontinuity will be specified later.

u = ū on S u , [19a] ũ = 0 on S u ,
Following standard procedures, the governing equations [18a] to [18c] can be cast in a weak form. To this end, the equilibrium equations [18a] are multiplied by the weight function w u , which is decomposed into ŵu and wu consistent with the displacement decomposition in [6], and integrated over the domain V to obtain a weak equilibrium statement which can be expanded using integration by parts, Gauss ' theorem and [18b]. More details on the derivation of the discrete weak governing equations can be found in [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] and [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF].

Discretized governing equations

Following a Bubnov-Galerkin approach, the linearized form of the discretized weak governing equation is obtained by substituting in the discrete weak governing equations the stress rate in the bulk

σ = D e ε = D e B ȧ + HB ḃ , [20]
with D e the tangent matrix for the bulk material, and the traction rate at a discontinuity

ṫ = T u = T N ḃ ⏐ ⏐ ⏐ S d , [ 21
]
where T relates traction rate ṫ and displacement jump rate u .

After standard manipulations, the linearized weak form of the governing equations reads (cf e.g. [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] and [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF])

K aa K ab K ba K bb δa δb = f ext,a f ext,b - f int,a f int,b , [22]
where

K aa = V B T D e B dV, [23a]
K ab = V + B T D e B dV, [23b] K ba = K T ab = V + B T D e B dV, [23c]
K bb = V + B T D e B dV + S d N T T N dS, [23d]
and

f int,a = V B T σ dV, f ext,a = S t N T t dS, [24a] f int,b = V + B T σ dV + S d N T t dS, f ext,b = S + t N T t dS. [24b]
In the above relations, B is a matrix containing derivatives of shape functions, σ is an array containing the components of the stress tensor in engineering notation and t is the array of the traction forces across the discontinuity.

Note that symmetry of the stiffness matrix is assured if the material tangent matrices D e and T are symmetric. As a final remark, the arrays N multiplying a and b in [9] are normally not the same since only part of the extra degrees of freedom in the array b might be activated. Consequently, the matrices B multiplying a and b in [20] are also normally not the same. However, since the system of equations in [22] is assembled only for the active degrees of freedom, it is still possible to use the standard N and B matrices. Care should of course be taken in considering only active degree of freedom in the strain computation. More details on the finite-element implementation can be found in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF], [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF][START_REF] Sukumar | Modeling quasi-static crack growth with the extended finite element method -Part I: Computer implementation[END_REF].

Finite element implementation

PU-based discontinuous elements require some important modifications to existing finite element codes. Some of the issues are illustrated in the example reported in Section 8. Here we make reference to some issues relevant to static and propagating discontinuities.

Numerical integration:

The integration of element matrices on part of the element domain requires the definition of special integration schemes. A common strategy consists in considering quadrature subcells aligned with the discontinuity line [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. When an element is intersected by a discontinuity, the two resulting do-discontinuity Figure 4. Composite integration scheme for a quadrilateral element crossed by a discontinuity (the crossed circles are integration points on the discontinuity and the black dots are integration points for the continuum) mains are triangulated and each triangular sub-domain is mapped to a parent unit triangle over which a three-point symmetric quadrature rule, with interior points within the triangular sub-domain, is considered (see Figure 4). The discontinuity contribution to the stiffness matrix and the internal force vector are integrated on the discontinuity through a two-point Gauss integration scheme. For elements not crossed by a discontinuity, a standard Gauss integration scheme is used.

It is worth noting that the resulting composite integration scheme for the continuum reported in Figure 4, which consists of 24 integration points in eight triangular sub-domains, integrates correctly a quadratic function over a parent quadrilateral element. In the numerical integration of the element matrices, this composite integration scheme results, for the terms containing the strain energy density, in a full integration rule for bilinear quadrilateral elements and in a reduced integration rule for a quadratic quadrilateral elements. The actual construction of the composite integration rule is based on an exact algebraic inverse iso-parametric mapping for bilinear quadrilateral elements. This technique can be used for quadratic quadrilateral elements provided that the element has straight sides and evenly spaced side nodes.

An alternative to the above procedure is to use a high order integration scheme in the elements so that there is at least one integration scheme in each domain. More recently, [START_REF] Ventura | On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method[END_REF] has proposed an approach in which discontinuous shape functions are replaced by equivalent polynomials thus eliminating the need for integration subcells.

Propagating discontinuities: Most commonly used criteria for discontinuity extension are based on stress intensity factors [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF], principal stresses [START_REF] Jirásek | Embedded crack model. Part II: Combination with smeared cracks[END_REF][START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] and loss of hyperbolicity [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF]. Criteria based on damage or plastic strain accumulation can be used [START_REF] Simone | From continuous to discontinuous failure in a gradientenhanced continuum damage model[END_REF][START_REF] Simone | The use of displacement discontinuities in a rate-dependent medium[END_REF][START_REF] Comi | An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation[END_REF] when the bulk material is inelastic.

As for the crack tip location, the two possible alternatives are to locate a crack tip within an element [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] or at its boundary [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF]. The latter approach has the advantage of greatly simplifying the computational implementation by setting the displacement jump at the discontinuity to zero. This is achieved by considering only standard a dofs for the nodes on an element boundary touched by a discontinuity (see Figure 5). When a discontinuity is extended in a neighboring element, the nodes behind the discontinuity tip are enhanced.

The direction of discontinuity extension is computed using the principal directions of a nonlocal stress tensor which is calculated as a weighted average of stresses using a Gaussian weight function [START_REF] Jirásek | Embedded crack model. Part II: Combination with smeared cracks[END_REF][START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] with an the interaction radius equal to three times the average element size ahead the discontinuity tip. The discontinuity propagates in the direction normal to the direction of the maximum nonlocal principal stress. When the discontinuity is close to a boundary, the discontinuity extension direction is aligned with the previous discontinuity segment to avoid an incorrect direction determination due to the bias produced by an unsymmetric domain in the weighting procedure.

Finally, to preserve the robustness of the Newton-Raphson solution procedure, a discontinuity is introduced as a straight segment at the end of a load increment in the element ahead of a discontinuity tip if the initiation criterion is fulfilled. This procedure is repeated in the elements ahead of this element with the extended discontinuity until the initiation criterion is no longer satisfied.

Many important details have been glossed over for the sake of brevity in this short summary; the interested reader is referred to [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], [START_REF] Dolbow | Discontinuous enrichment in finite elements with a partition of unity method[END_REF], [START_REF] Sukumar | Modeling quasi-static crack growth with the extended finite element method -Part I: Computer implementation[END_REF], and [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF] for further details.

A minimal one-dimensional example

Some of the issues related to the implementation of PU-based discontinuous elements are illustrated by means of a simple one-dimensional example of a bar in tendiscontinuity tip sion as depicted in Figure 6. The solution will be given for standard and for PU-based discontinuous elements. The bar is modelled by means of the two discretizations depicted in Figure 7. The spring is first represented by a conventional interface element (Figure 7a) while the discontinuous interpolation of [9] is exploited for the discretization depicted in Figure 7b. The domain V + is defined by 0 < x < L in Figure 6.

In the following, E is the Young's modulus, A is the cross-sectional area and d is the spring stiffness. The element types needed for the finite-element solution are depicted in Figure 8.

Stiffness matrix computation.

The stiffness matrix of the one-dimensional truss element of length L depicted in Figure 9a is

K truss = EA L 1 -1 -1 1 .
[25]

The same stiffness matrix describes a one-dimensional conventional interface element which can be conceived as a translational spring element with stiffness d = EA/L in [25].

The sub-matrices in [23] are expanded for the truss element with a discontinuity in the middle section depicted in Figure 9b. Note that the positive part V + of the domain goes from x = 0 to x = L/2 (H = 1 for x < L/2; see Figure 9b) and the presence of extra degrees of freedom b for both nodes. Sub-matrix K aa is the same as K truss . The remaining sub-matrices are expanded as In K bb , the volume integral equals K ab while the surface integral is evaluated on x = L/2 and yields

K ab = K ba = V + B T D e B dV = L/2 0 B T D e B dx = EA 2L 1 -1 -1 1 . [26] 1 3 4 L L 1 L L/2 2 3 4 L/2 (
1 2 4 3 1 4 2 3 (a) (b) 2 3 (c)
S d N T T N dS = N T T N ⏐ ⏐ ⏐ x=L/2 = d 4 1 1 1 1 . [27]
Assembly of the sub-matrices into the element stiffness matrix yields If the discontinuity does not cross the element but extra degrees of freedom b are active at one of the element nodes (see e.g. Figure 8a), the element stiffness matrix can be derived following the procedure described above and reads as

K = K aa K ab K ba K bb = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ EA L -EA L EA 2L -EA 2L -EA L EA L -EA 2L EA 2L EA 2L -EA 2L EA 2L + d 4 -EA 2L + d 4 -EA 2L EA 2L -EA 2L + d 4 EA 2L + d 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . [28] j L i a i a j x i j L x a i a j (a) (b) b i b j V - V + S d
K = 2EA L ⎡ ⎣ 1 -1 -1 -1 1 1 -1 1 1 ⎤ ⎦ , [29]
where the degrees of freedom have been ordered in the sequence a 1 a 2 b 2 (note that the truss length is L/2).

Assembly and solution.

For the bar with the spring modelled by a conventional interface element (see Figure 7a), assembly of local stiffness matrices into the global stiffness matrix for the active degrees of freedom results in the system of equations For the discretization with the PU-based discontinuous element (see Figure 7b), the global system of equations reads

⎡ ⎢ ⎢ ⎣ EA L + d -d 0 -d EA L + d -EA L 0 -EA L EA L ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ a 2 a 3 a 4 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ 0 0 P ⎤ ⎥ ⎥ ⎦ , [30]
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 3EA L -EA L 5 2 EA L -1 2 EA L 0 -EA L 3EA L -1 2 EA L 1 2 EA L -2EA L 5 2 EA L -1 2 EA L 5 2 EA L + d 4 -1 2 EA L + d 4 0 -1 2 EA L 1 2 EA L -1 2 EA L + d 4 1 2 EA L + d 4 0 0 -2EA L 0 0 2EA L ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a 2 a 3 b 2 b 3 a 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 P ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
and its solution yields For interface elements, the displacement jump across the discontinuity S d is expressed as the difference u h = a 3a 2 = P/d of the displacement of the doubled nodes. In PU-based discontinuous elements, the displacement jump is expressed through [10] as u h = N b| S d = 0.5 × (b 2 + b 3 ) = -P/d (see Figure 10). Note that the minus sign in the displacement jump value indicates that the local frame of the discontinuity is in the opposite direction with respect to the global frame x-in other words, the negative sign is due to the choice of the positive part of the domain. Figure 11. Occurrence of oscillations in the traction profile for some values of the interface stiffness and with Gauss integration scheme for the interface contribution (from Simone ( 2004))

PU-based discontinuous elements and the "dummy stiffness" approach

The "dummy stiffness" approach is often used in combination with conventional interface elements to simulate perfect contact prior to the loss of material coherence along potential crack lines (like e.g. the dotted line in Figure 11a). This approach consists in considering a high value of the interface stiffness and is known to produce unrealistic oscillations in the traction profile like those shown in Figure 11b [START_REF] Schellekens | On the numerical integration of interface elements[END_REF]. These oscillations can however be eliminated by an appropriate choice of the interface integration scheme.

PU-based discontinuous elements, when used with a "dummy stiffness" approach, suffer from the same type of problem [START_REF] Simone | Partition of unity-based discontinuous elements for interface phenomena: Computational issues[END_REF]. Similar to classical interface elements, oscillations in the traction profile are linked to the appearance of pathological coupling between degrees of freedom in the part of the stiffness matrix responsible for the "interface" behavior (matrix K bb in [23]b). Unlike interface elements, these oscillations are not removed by using a particular integration scheme; they can be however circumvented by activating the degrees of freedom responsible for the displacement jump when they are required. This problem may still be present if favorable conditions exist, e.g. high stress gradient and high discontinuity stiffness in elastic regime or cyclic loading conditions which promote crack closure. In these cases, it is worth considering the solution proposed by [START_REF] Mourad | A bubble-stabilized finite element method for Dirichlet constraints on embedded interfaces[END_REF].

Crack propagation in a single-edge notched beam

A single-edge notched beam, depicted in Figure 12, is subjected to an antisymmetric four-point-shear loading [START_REF] Schlangen | Experimental and Numerical Analysis of Fracture Processes in Concrete[END_REF] right support. The beam is in concrete with Young's modulus E = 35 GPa, Poisson's ratio ν = 0.2, tensile strength f t = 3.0 MPa and fracture energy G f = 0.1 N/mm. The material parameters have been taken from [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF]. The loading platens have a Young's modulus one order of magnitude larger than the concrete. The load is applied by means of an indirect displacement control procedure with the crack mouth sliding displacement (cmsd), defined as the relative vertical displacement of the opposite faces of the notch, taken as control parameter. The boundary conditions are specified by constraining the displacement at the upper-right support in both directions and by constraining the vertical displacement at the upper-left support. The beam is analyzed under plane strain conditions. The simulation is performed using 1184 six node triangular elements. Crack propagation has been obtained using a principal stress initiation criterion. In this application the bulk material is considered as linear elastic whereas the normal traction-separation law at the discontinuity is defined by

t n = f t exp - f t G f κ , [33]
where f t is the tensile strength of the material and G f is the fracture energy. κ is a history parameter equal to the largest value of the normal separation at the interface. Shear traction contributions can be incorporated in the formulation if necessary. In this example [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] used

f 2 t G f exp -f t G f κ 0 h s d int exp (h s κ) u s d int exp (h s κ) u n u s . [35]
The predicted crack path shown in Figures 13 and 14 is curved and is in agreement with the experimental results [START_REF] Schlangen | Experimental and Numerical Analysis of Fracture Processes in Concrete[END_REF]. Worth noting is that the crack is independent of the mesh (see close-up of crack path in Figure 13).

In terms of global quantities, the numerical load-cmsd curves are compared to the experimental response in Figure 15. The agreement is quite good for both constant and variable shear stiffness up to the first part of the post-peak branch; the numerical curves are slightly too brittle for cmsd > 0.3 mm indicating that a more sophisticated expression of the traction-separation law is needed. [START_REF] Schlangen | Experimental and Numerical Analysis of Fracture Processes in Concrete[END_REF] versus numerical [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] load-crack mouth sliding displacement (cmsd) for variable shear stiffness (d int = 1.0 × 10 2 N/mm, d κ=1.0 = 1.0 × 10 -4 N/mm) and constant shear stiffness (d int = 1.0 N/mm)

Concluding remarks

Partition of unity-based discontinuous elements are a powerful alternative to conventional interface elements in problems with evolving boundaries. One of the obvious advantages, compared to conventional interface elements, is that mesh adjustment along the propagating boundary is avoided since a discontinuity can pass arbitrarily through solid elements. This extremely flexible and powerful approach has been applied to a wide variety of problems which could not be described here. For more information on the partition of unity method and its extension to discontinuous element representation, the interested reader is referred to the publications that reference Duarte et al. (Duarte et al., 1996b;[START_REF] Duarte | Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems[END_REF]Duarte et al., 1996a;[START_REF] Oden | A new cloud-based hp finite element method[END_REF], [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF], and [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

As a final remark, it bears emphasis that XFEM, PUFEM and GFEM are examples of the so-called partition of unity method which originated in the works of Babuška et al. (Babuška et al., 1994;[START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] and Duarte and Oden [START_REF] Duarte | Hp Clouds-A meshless method to solve boundary-value problems[END_REF]Duarte et al., 1996b).

NOTE. -This paper is mainly based on Duarte et al. (2007) and [START_REF] Simone | Partition of unity-based discontinuous elements for interface phenomena: Computational issues[END_REF].

Alfano G., Crisfield M. A., « Finite element interface models for the delamination analysis of laminated composites: Mechanical and computational issues », International Journal for Numerical Methods in Engineering, vol. 50, n°7, 2001Engineering, vol. 50, n°7, , p. 1701Engineering, vol. 50, n°7, -1736. . Babuška I., Caloz G., Osborn J. E., « Special finite element methods for a class of second order elliptic problems with rough coefficients », SIAM Journal on Numerical Analysis, vol. 31, n°4, 1994, p. 945-981.

Figure 1 .

 1 Figure 1. A GFEM shape function φ Ii as the product of a partition of unity N I and an enrichment function E i (Duarte et al., 2006) (the enrichment takes the form of a polynomial function in (a) and of a discontinuous function in (b))

Figure 3 .

 3 Figure 3. Nodes affected by enrichment in a crack problem

Figure 5 .Figure 6 .

 56 Figure 5. Extension of a discontinuity within the body V (the black circles indicates nodes with standard degrees of freedom a while the white circles indicates nodes with extra degrees of freedom b for the nodes whose support is crossed by a discontinuity)

Figure 7 .

 7 Figure 7. Discretization for the tension test with (a) conventional interface elements and (b) PU-based discontinuous elements (the black circles indicates the standard degrees of freedom a while the white circles indicates the extra degrees of freedom b; the discontinuity is indicated by the vertical dotted line).

Figure 8 .

 8 Figure 8. Element assembly for the discretization with PU-based discontinuous elements (note the effect of the choice of the domain V + on the activation of the extra degrees of freedom b)

Figure 9 .

 9 Figure 9. One-dimensional (a) truss element and (b) PU-based discontinuous truss element

Figure 10 .

 10 Figure10. Displacement field for the one-dimensional tension test: (a) total displacement field given by the sum of (b) the û field and (c) the Hũ field (see[32] for the numerical values)

Figure 12 .

 12 Figure 12. Single-edge notched beam (depth = 100 mm; all dimensions in mm).

Figure 13 .

 13 Figure 13. Curved crack path in the SEN beam (from Wells et al. (2001))

Figure 14 .

 14 Figure 14. Deformation of the SEN beam and σ xx stresses. The displacements are amplified by a factor 100 (from Remmers (2006))

t

  s = d int exp (h s κ) u s[34] with h s = ln d κ=1.0 d int , d int the initial crack shear stiffness, u s the crack sliding displacement, and d κ=1.0 the crack shear stiffness when κ = 1.0. The material tangent T ( ṫ = T u ) can be expressed by differentiating the expressions of the normal and tangential traction:

Figure 15 .

 15 Figure 15. Experimental[START_REF] Schlangen | Experimental and Numerical Analysis of Fracture Processes in Concrete[END_REF] versus numerical[START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF] load-crack mouth sliding displacement (cmsd) for variable shear stiffness (d int = 1.0 × 10 2 N/mm, d κ=1.0 = 1.0 × 10 -4 N/mm) and constant shear stiffness (d int = 1.0 N/mm)