# An eXtended Stochastic Finite Element Method for solving stochastic partial differential equations on random domains

Anthony Nouy, Alexandre Clément, Franck Schoefs, Nicolas Moes

► To cite this version:

Anthony Nouy, Alexandre Clément, Franck Schoefs, Nicolas Moes. An eXtended Stochastic Finite Element Method for solving stochastic partial differential equations on random domains. Computer Methods in Applied Mechanics and Engineering, 2008, 197 (51-52), pp.4663-4682. 10.1016/j.cma.2008.06.010. hal-00366617

HAL Id: hal-00366617

https://hal.science/hal-00366617v1

Submitted on 9 Mar 2009

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An eXtended Stochastic Finite Element Method for solving stochastic partial differential equations on random domains

A. Nouy*, A. Clément, F. Schoefs, N. Moës

Research Institute in Civil Engineering and Mechanics (GeM), Nantes Atlantic University, Ecole Centrale Nantes, UMR CNRS 6183, 2 rue de la Houssinière, B.P. 92208, 44322 Nantes Cedex 3, FRANCE

## Abstract

Recently, a new strategy was proposed to solve stochastic partial differential equations on random domains. It is based on the extension to the stochastic framework of the eXtended Finite Element Method (X-FEM). This method leads by a “direct” calculus to an explicit solution in terms of the variables describing the randomness on the geometry. It relies on two major points: the implicit representation of complex geometries using random level-set functions and the use of a Galerkin approximation at both stochastic and deterministic levels. In this article, we detail the basis of this technique, from theoretical and technical points of view. Several numerical examples illustrate the efficiency of this method and compare it to other approaches.

**Key words:** Computational Stochastic Mechanics, Stochastic Partial Differential Equations, Random domain, eXtended Finite Element Method, Stochastic Finite Element, Polynomial Chaos

## 1 Introduction

Computer simulations of mechanical models, supported by the availability of enhanced computational resources, take today a very significant place in decision-making which can have major consequences in economic or human terms. In a structural analysis, the incorporation of uncertainties inherent in

<footnote>* Corresponding author. Tel.: +33(0)2-51-12-55-20; Fax: +33(0)2-51-12-52-52
Email address: anthony.nouy@univ-nantes.fr (A. Nouy).
URL: http://www.univ-nantes.fr/nouy-a (A. Nouy).</footnote>
the model, related to material properties, loadings or geometry, seems today
essential if one seeks to obtain “reliable” numerical predictions, usable in a
design process or a decision-making. This necessity led to a rapid development
of many ad hoc numerical methods, such as stochastic finite elements methods
[1–3]. These methods provide high quality predictions which are explicit in
terms of the random variables describing the uncertainties. The increase in
their performances, due to the development of ad hoc resolution techniques
[4,6–8], allows expecting their use in many domains of applications such as
structural design, reliability analysis, etc. The incorporation of uncertainties
on material properties or loadings is quite well mastered within the framework
of these techniques. However, there is still no available efficient strategy to deal
with uncertainties on the geometry although it could have a great interest in
various applications.

A natural way to solve a stochastic problem defined on a random domain con-
sists in using a classical stochastic finite element method with remeshings, such
as Monte-Carlo simulations [9,10], response surface method, projection [11] or
regression [12] methods. These techniques only require the use of a simple
deterministic finite element calculation code. However, they require numerous
deterministic computations, each of which requiring the construction of a new
conforming mesh, which leads to prohibitive computational costs.

In order to avoid remeshings and to obtain an “explicit” description of the
stochastic solution, a possible alternative technique consists in building a clas-
sical finite element approximation on a reference deterministic domain, and to
introduce a random mapping between this domain and the random domain.
This technique can also be seen as a classical finite element technique on a
random mesh. Such a strategy has been used in [13,14] and also in [5] for
the case of random interfaces in layered media. The main difficulty lies in the
construction of a suitable random mapping in the case of complex geometries.

Recently, a new stochastic finite element method was proposed for solving
partial differential equations on random domains [15]. This technique, called
eXtended Stochastic Finite Element method (X-SFEM), is based on an ex-
tension to the stochastic framework of the X-FEM method [16,17]. It relies
on two major points. The first one consists in representing the geometry in
an implicit way by using the level-set technique [18]. The random domain is
then characterized by the negative set of a random level-set function, which
is defined on a fictitious deterministic domain containing all outcomes of the
physical domain. The second point consists in using a Galerkin approxima-
tion technique at both deterministic and stochastic levels. A tensor product
approximation is made possible by considering prolongation of functions on
the fictitious domain. The resulting technique does not require any remeshing
and allows the handling of complex geometries. Another advantage is that it
leads to a solution which is explicit in terms of the basic random variables de-
scribing the geometry, thus allowing an easy post-processing of the stochastic
solution at a very low cost.

In this paper, we introduce the basis of the X-SFEM method from a theo-
retical point of view and also focus on technical aspects. The method will be
presented in the context of linear elasticity. However, the mathematical results
and the proposed methodology could be applied to more general linear ellip-
tic stochastic partial differential equations defined on random domains. Here,
we only deal with random shapes. This is a degenerated case of the X-FEM
method [19,20] where only a classical finite element approximation can be
used at the space level, without enrichment of the approximation space by the
partition of unity method [21]. Let us mention that fictitious domain methods
[22–25], recently extended to the stochastic framework [26], start with similar
reformulation of the problem and definition of approximation spaces. How-
ever, the methodology proposed in this article differs by the representation
of the geometry and the underlying computational aspects (construction and
resolution of the discretized problem).

The outline of the paper is as follows. In section 2, we introduce the strong and weak formulations of a linear elasticity problem defined on a random domain. In section 3, we briefly recall the basis of stochastic finite element methods in the case of a deterministic domain. It allows us to point out why their extension in the case of a random domain is not straightforward and also to introduce some useful notations. In section 4, we present possible finite element strategies to deal with random domains and point out some of their drawbacks. Section 5 is devoted to the presentation of the X-SFEM method from a theoretical point of view. In particular, we focus on the mathemati-cal interpretation of the approximation based on a Galerkin projection and compare it to other projection methods. Section 6 details the computational aspects of the X-SFEM method: construction and resolution of the discretized problem. In particular, we illustrate that the construction of the discretized problem requires the development of a specific numerical integration technique at the stochastic level. For that purpose, a general numerical procedure is pro-posed, which consists in building suitable partitions of the stochastic domain and using a Gaussian quadrature on each subdomain. Finally, in section 7, three numerical examples illustrate several aspects of the X-SFEM method: quality of the Galerkin approximation, influence of the proposed stochastic integration procedure and efficiency of the resolution strategy. The method is systematically compared to other approaches. In particular, it is compared to a stochastic L² projection method based on X-FEM.
## 2 Formulation of the problem

### 2.1 The deterministic problem

In this article, we will focus on the analysis of the deformation of a structure under small perturbations assumption. We consider that the structure occupies a domain $\Omega \subset \mathbb{R}^d$ (see figure 1). It is submitted to body loads $f(x)$ on $\Omega$ and surface loads $F(x)$ on a part $\Gamma_2$ of the boundary $\partial\Omega$. The structure is fixed on a part $\Gamma_1$ of the boundary $\partial\Omega$ such that $\Gamma_1 \cap \Gamma_2 = \emptyset$. The complementary part of $\Gamma_1 \cup \Gamma_2$ in $\partial\Omega$, denoted by $\Gamma_0$, is considered as a free boundary. We denote by $u(x)$ the displacement field, $\varepsilon(u)$ the symmetric part of the displacement gradient (i.e. the strain tensor) and by $\sigma$ the stress tensor. In this article, we will consider the case of a linear elastic material and denote by $C$ the Hooke fourth-order tensor. The strong formulation of the problem writes: find $(u, \sigma)$ such that

$$ 
\begin{aligned}
& \mathit{div} \, \sigma + f = 0 \quad \text{on } \Omega, \\
& \sigma = C : \varepsilon(u) \quad \text{on } \Omega, \\
& \sigma \cdot n = F \quad \text{on } \Gamma_2, \\
& \sigma \cdot n = 0 \quad \text{on } \Gamma_0, \\
& u = 0 \quad \text{on } \Gamma_1.
\end{aligned}
\tag{1} $$

where $n$ is the unit outward normal to the boundary.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Model problem</figcaption></figure>

### 2.2 The stochastic problem

Now, we consider that the structure occupies a random domain. For the sake of generality, we also consider that material properties and loads can be random. Let us suppose that the probability space $(\Theta, \mathcal{B}, P)$ allows the representation
of the probabilistic content of the problem. $\Theta$ is the set of elementary events, $\mathcal{B}$ a $\sigma$-algebra on $\Theta$ and $P$ a probability measure. An element $\theta \in \Theta$ is an elementary event. The random domain can be characterized by the random variable $\Omega: \theta \in \Theta \mapsto \Omega(\theta) \subset \mathbb{R}^d$. $\Omega(\theta)$ represents an outcome of the domain. Figure 2 illustrates two outcomes of the random domain and of the associated parts of its boundary.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Two outcomes of the random domain $\Omega(\theta)$ and associated parts of boundary $\Gamma_1(\theta)$ and $\Gamma_2(\theta)$.</figcaption></figure>

Then, the strong formulation of the stochastic problem writes: find random fields $u$ and $\sigma$ which verify $P$-almost surely the set of equations (1). We recall that a property is $P$-almost sure if the probability of the event “the property is false” is zero.

## 2.3 Weak formulation of the stochastic problem

Let us first introduce a classical weak formulation at the space level: find $u$ such that $\forall \theta \in \Theta$, $u(\cdot, \theta) \in \mathcal{U}(\theta) = \{v \in (H^1(\Omega(\theta)))^d; v = 0 \text{ on } \Gamma_1(\theta)\}$, and such that we have $P$-almost surely

$$a(u(\cdot, \theta), v; \theta) = b(v; \theta) \quad \forall v \in \mathcal{U}(\theta), \qquad (2)$$

where $a(\cdot, \cdot; \theta)$ is a continuous symmetric coercive bilinear form on $\mathcal{U}(\theta)$ and $b(\cdot; \theta)$ is a continuous linear form on $\mathcal{U}(\theta)$, defined by

$$a(u, v; \theta) = \int_{\Omega(\theta)} \varepsilon(v) : C : \varepsilon(u) \, dx, \qquad (3)$$

$$b(v; \theta) = \int_{\Omega(\theta)} f \cdot v \, dx + \int_{\Gamma_2(\theta)} F \cdot v \, ds. \qquad (4)$$

Then, a weak formulation at both space and stochastic levels can be introduced. We consider that the problem can be formulated in the following func-
tion space:

$$
\begin{align}
\mathcal{V} &= \{ \mathbf{u} : \theta \in \Theta \mapsto \mathbf{u}(\cdot, \theta) \in \mathcal{U}(\theta); \int_{\Theta} \|\mathbf{u}\|_{\mathcal{U}(\theta)}^2 dP(\theta) < \infty \} \nonumber \\
&:= L^2(\Theta, dP; \mathcal{U}). \tag{5}
\end{align}
$$

The weak formulation then writes : find **u** ∈ V such that

$$
A(\mathbf{u}, \mathbf{v}) = L(\mathbf{v}) \quad \forall \mathbf{v} \in \mathcal{V}, \tag{6}
$$

where

$$
A(\mathbf{u}, \mathbf{v}) = \int_{\Theta} a(\mathbf{u}, \mathbf{v}; \theta) dP(\theta) := E(a(\mathbf{u}, \mathbf{v}; \theta)), \quad (7)
$$

$$
L(\mathbf{v}) = \int_{\Theta} b(\mathbf{v}; \theta) dP(\theta) := E(b(\mathbf{v}; \theta)). \quad (8)
$$

$E(\cdot)$ denotes the mathematical expectation. We will consider that bilinear
form $A$ is continuous and coercive on $\mathcal{V}$ and that linear form $L$ is continuous
on $\mathcal{V}$, such that problem (6) is well posed. In the case of a deterministic
domain, the reader can refer to [2,27] for a detailed and comprehensive study
of the regularity requirements on the data ($C$, $f$ and $F$) leading to the desired
properties. To the knowledge of the authors, stochastic problem (6) has not
been deeply studied from a mathematical point of view in the case of random
domains. In particular, regularity requirements on the random domain should
be further analyzed. These mathematical aspects are beyond the scope of this
article. Some mathematical results on the effects of uncertainties in the domain
can be found in [28,29].

**Remark 1** Coercivity and continuity properties of A and L are inherited from those of a and b: ∀**u**, **v** ∈ U(θ), |a(**u**, **v**; θ)| ≤ C_a||**u**||_{U(θ)}||**v**||_{U(θ)}, a(**v**, **v**; θ) ≥ α_a||**v**||_{U(θ)}^2, and |b(**v**; θ)| ≤ C_b||**v**||_{U(θ)}, where C_a, C_b < ∞ and α_a > 0. These three constants depends on material properties, loads and also on the geometry. If there exists such constants C_a, C_b and α_a independent of the event θ, coercivity and continuity properties of A and L follows with the same constants.

2.4 Stochastic modeling and discretization

In practice, we generally consider that the probabilistic content of the stochas-
tic problem can be represented by a finite set of random variables ξ : Θ → Rm.
This is the case when material properties, loads and geometry depend on a
finite set of parameters which are random variables. When these parameters
are stochastic fields, a preliminary step is usually performed in order to dis-
cretize these parameters (e.g. by Karhunen-Loeve decomposition [30,31,1])
and represent them in terms of a finite set of random variables.
Let us now consider that the data of stochastic problem (6) can be written in terms of $\xi$, so that bilinear form $a$ and linear form $b$ can also be written in terms of $\xi$: $a(\cdot, \cdot; \theta) \equiv a(\cdot, \cdot; \xi(\theta))$, $b(\cdot; \theta) \equiv b(\cdot; \xi(\theta))$. Then, by the Doob-Dynkin's lemma [32], we have that the solution of problem (6) can be written in terms of $\xi$, i.e. $u(\cdot, \theta) \equiv u(\cdot, \xi(\theta))$. The stochastic problem can then be reformulated in the $m$-dimensional image probability space $(\Theta, \mathcal{B}, P_{\xi})$, where $\Theta \subset \mathbb{R}^m$ is the range of $\xi$ and $P_{\xi}$ the probability measure associated with $\xi$ (see e.g. [2] for details). Function space (6) is now replaced by

$$ \mathcal{V} = \{u : \xi \in \Theta \mapsto u(\cdot, \xi) \in \mathcal{U}(\xi); E(\|u(\cdot, \xi)\|_{\mathcal{U}(\xi)}^2) < \infty\}. \quad (9) $$

In the definition of function space $\mathcal{V}$ and of bilinear form $A$ and linear form $L$, the mathematical expectation must be interpreted as $E(f(\xi)) = \int_{\Theta} f(y)dP_{\xi}(y)$.

In what follows, we will mainly use the initial probability space $(\Theta, \mathcal{B}, P)$. The reader must keep in mind that at each moment, the elementary event $\theta \in \Theta$ can be replaced by $\xi \in \Theta$ in the expression of any random function.

### 3 Stochastic finite element methods in the case of a deterministic domain

Before introducing some classical techniques to deal with random domains, we briefly recall the basis of stochastic finite element methods in the case of a deterministic domain (see e.g. [1–3,33]). It will allow us to introduce some useful notations and to point out why their extension in the case of a random domain is not straightforward.

Here, we suppose that the domain $\Omega$ and the associated parts of its boundary $\Gamma_i$ are deterministic. In this case, function space $\mathcal{U}$ does not depend on the elementary event $\theta$. Then, function space $\mathcal{V}$ has a tensor product structure:

$$ \mathcal{V} = L^2(\Theta, dP; \mathcal{U}) \cong L^2(\Theta, dP) \otimes \mathcal{U} := \mathcal{S} \otimes \mathcal{U}, \quad (10) $$

where $\mathcal{S} = L^2(\Theta, dP)$ is the space of second order real-valued random variables. Classical stochastic finite elements methods then introduce a tensor product approximation as follows.
## 3.1 Approximation at the space level

At the space level, we introduce a finite element approximation space $\mathcal{U}_h$ defined by

$$ \mathcal{U}_h = \{ \mathbf{v}(\mathbf{x}) = \sum_{i=1}^{N} \varphi_i(\mathbf{x}) v_i = \varphi(\mathbf{x}) \mathbf{v}, \mathbf{v} \in \mathbb{R}^N \}, \quad (11) $$

where the $\varphi_i(\mathbf{x}) \in \mathcal{U}$ are the finite element basis functions associated with a mesh $\mathcal{T}_h$ of domain $\Omega$. A semi-discretized approximation space $\mathcal{V}_h \subset \mathcal{V}$ can then be defined as follows:

$$ \begin{aligned} \mathcal{V}_h &= \mathcal{S} \otimes \mathcal{U}_h \\ &= \{ \mathbf{v}(\mathbf{x}, \theta) = \sum_{i=1}^{N} \varphi_i(\mathbf{x}) v_i(\theta) = \varphi(\mathbf{x}) \mathbf{v}(\theta), \mathbf{v} \in \mathcal{S} \otimes \mathbb{R}^N \}. \end{aligned} \quad (12) $$

The semi-discretized approximate solution $\mathbf{u}_h \in \mathcal{V}_h$ of problem (6) is defined by:

$$ A(\mathbf{u}_h, \mathbf{v}_h) = L(\mathbf{v}_h) \quad \forall \mathbf{v}_h \in \mathcal{V}_h. \quad (13) $$

Problem (13) can be written under the following matrix strong form: find $\mathbf{u}: \theta \mapsto \mathbf{u}(\theta) \in \mathbb{R}^N$ such that we have P-almost surely

$$ \mathbf{A}(\theta)\mathbf{u}(\theta) = \mathbf{b}(\theta), \quad (14) $$

where random matrix $\mathbf{A}$ and random vector $\mathbf{b}$ are defined by:

$$ (\mathbf{A}(\theta))_{ij} = a(\varphi_j, \varphi_i; \theta), \quad (15) $$

$$ (\mathbf{b}(\theta))_i = b(\varphi_i; \theta). \quad (16) $$

## 3.2 Approximation at the stochastic level

Several choices have been proposed for the construction of a basis of the space of second order random variables $\mathcal{S} = L^2(\Theta, dP)$: polynomial chaos [1,34], generalized polynomial chaos [35,36], h-p finite elements [37] or wavelets [38,39]. Let $\{H_\alpha\}_{\alpha \in \mathcal{J}}$ be such a basis, where $\mathcal{J}$ denotes a countable set of indices. Then, we can naturally introduce a finite dimensional approximation space $\mathcal{S}_P \subset \mathcal{S}$, defined by

$$ \mathcal{S}_P = \{ v(\theta) = \sum_{\alpha \in \mathcal{J}_P} v_\alpha H_\alpha(\theta), v_\alpha \in \mathbb{R} \}, \quad (17) $$

where $\mathcal{J}_P$ is a subset of $P$ indices in $\mathcal{J}$. In the following, we will suppose that the $H_\alpha$ are orthonormal, i.e. $E(H_\alpha H_\beta) = \delta_{\alpha\beta}$.
Let us illustrate the construction of the generalized polynomial chaos approximation [35,36], which will be used in numerical examples. We consider, as in section 2.4, that the probabilistic content of the problem is represented by m random variables ($ξ_1, ..., ξ_m$) := $ξ$. For the construction of chaos basis associated with an arbitrary probability measure, the reader can refer to [36]. Here, we will consider that random variables are statistically independent (property possibly obtained by a suitable mapping). Let ($\Theta$, $\mathcal{B}$, $P_ξ$) be the associated image probability space. The space $L^2(\Theta, dP_ξ)$ of second order random variables on $\Theta$ has a tensor product structure: $L^2(\Theta; dP_ξ) = \bigotimes_{j=1}^m L^2(\Theta^{(j)}, dP_{\xi_j})$, where $\Theta^{(j)} \subset \mathbb{R}$ is the range of random variable $ξ_j$ and $P_{\xi_j}$ the associated marginal probability measure. Basis functions $H_α$ are then taken as multi-dimensional orthonormal polynomials. A simple way to define the $H_α$ is to use a tensor product basis. Let us denote by $\mathcal{J} = \{α = (α_1, ..., α_m) ∈ \mathbb{N}^m\}$ the set of m-dimensional multi-indices. Then, denoting by {$h_k^{(j)}$}$_{k≥0}$ the set of orthonormal polynomials in $L^2(\Theta^{(j)}, dP_{\xi_j})$, where index k denotes the degree of the polynomials, the $H_α$ are simply defined as follows: $H_α(ξ) = \prod_{j=1}^m h_{\alpha_j}^{(j)}(ξ_j)$. The approximation space $S_P$ can then be defined by selecting the multi-dimensional polynomials with total degree less than p. It consists in taking the following subset of indices $\mathcal{J}_P = \{α ∈ \mathcal{J}; \sum_{i=1}^m α_i ≤ p\}$. In this case, the number of basis functions is $P = \frac{(m+p)!}{m!p!}$.

Finally, a tensor product approximation space $V_{h,P} \subset \mathcal{V}$ is classically defined as follows:

$$ 
\begin{aligned}
V_{h,P} &= S_P \otimes U_h \qquad (18) \\
&= \{v(x, \theta) = \sum_{i=1}^{N} \sum_{\alpha \in \mathcal{J}_P} \varphi_i(x) H_\alpha(\theta) v_{i,\alpha}, v_{i,\alpha} \in \mathbb{R}\} \\
&= \{v(x, \theta) = \varphi(x) v(\theta), v(\theta) \in S_P \otimes \mathbb{R}^N\}.
\end{aligned}
$$

### 3.3 Definition of the approximation

Many choices are possible for the definition of an approximation in $V_{h,P}$. Below, we present the more commonly used definitions, namely the stochastic Galerkin method and the classical projection method.

#### 3.3.1 Stochastic Galerkin method

A classical way to define the approximate solution $u_{h,P} \in V_{h,P}$ consists in injecting the approximation space $V_{h,P}$ in the weak formulation (6). It leads
to the so-called Galerkin approximation which solves

$$A(\mathbf{u}_{h,P}, \mathbf{v}_{h,P}) = L(\mathbf{v}_{h,P}) \quad \forall \mathbf{v}_{h,P} \in \mathcal{V}_{h,P}, \qquad (19)$$

which is equivalent to the following system of $N \times P$ equations:

$$\sum_{\alpha \in \mathcal{J}_P} E(H_{\beta} \mathbf{A} H_{\alpha}) \mathbf{u}_{\alpha} = E(\mathbf{b} H_{\beta}) \quad \forall \beta \in \mathcal{J}_P. \qquad (20)$$

### 3.3.2 Projection method

The projection method consists in defining the approximate solution as the projection on $\mathcal{V}_{h,P}$ of the semi-discretized solution with respect to a suitable inner product. A classical choice consists in using the natural $L^2$ inner product on $\mathcal{V}$, defined by

$$\ll \mathbf{u}, \mathbf{v} \gg_{L^2} = \int_{\Theta} \int_{\Omega} \mathbf{u} \cdot \mathbf{v} \, dx \, dP(\theta). \qquad (21)$$

Denoting by $\mathbf{u}_h(\mathbf{x}, \theta) = \varphi(\mathbf{x})\mathbf{u}(\theta) \in \mathcal{V}_h$ the semi-discretized solution, which solves (13), the approximation $\mathbf{u}_{h,P} \in \mathcal{V}_{h,P}$ is then defined by:

$$\ll \mathbf{u}_{h,P}, \mathbf{v}_{h,P} \gg_{L^2} = \ll \mathbf{u}_h, \mathbf{v}_{h,P} \gg_{L^2} \quad \forall \mathbf{v}_{h,P} \in \mathcal{V}_{h,P}. \qquad (22)$$

Equation (22) is equivalent to the following definition of the coefficients of the approximation:

$$\mathbf{u}_{\alpha} = \int_{\Theta} H_{\alpha}(\theta) \mathbf{u}(\theta) dP(\theta) = E(H_{\alpha} \mathbf{u}). \qquad (23)$$

We emphasize on the fact that this equivalence comes from the tensor product structure of the approximation and of the particular choice of inner product. The expectation can then be computed by a suitable numerical integration at the stochastic level (quadrature, Monte-Carlo, etc.), which writes

$$\mathbf{u}_{\alpha} \approx \sum_{k} \omega_k H_{\alpha}(\theta_k) \mathbf{u}(\theta_k), \qquad (24)$$

where the $\theta_k$ are the integration points (particular elementary events) and the $\omega_k$ the associated integration weights. For each integration point, computing $\mathbf{u}(\theta_k) = \mathbf{A}^{-1}(\theta_k)\mathbf{b}(\theta_k)$ requires solving a classical deterministic problem. In this case, the projection method is said to be non-intrusive since it only requires the use of a simple deterministic code for each integration point. However, computational costs are generally prohibitive since this technique often necessitates a large number of integration points in order to obtain a good accuracy on the coefficients. Let us mention that sparse quadrature techniques have been proposed in order to reduce the number of integration points [3].
# 4 Classical stochastic finite elements to deal with random domains

The question is now: how to define an approximation of problem (6) in the case of a random domain? The extension of classical stochastic finite element methods presented in section 3 is not straightforward since $\mathcal{V}$, defined in (5), has no more a tensor product structure. In the following, we recall two possible strategies to construct stochastic finite element approximations in the case of a random domain.

## 4.1 Classical finite element techniques with remeshings

A natural way to solve a stochastic problem defined on a random domain consists in using a classical stochastic finite element method with remeshings, such as Monte-Carlo simulation [9,10], response surface method, projection [11] or regression [12] methods. These techniques require solving problem (2) for given elementary events $\theta_k \in \Theta$. The meaning of these points depends on the method which is used: random samplings for Monte-Carlo, stochastic integration points for projection method, etc. These techniques have the great advantage that they only require the use of a simple deterministic finite element calculation code.

For a given elementary event $\theta_k \in \Theta$, we introduce a new finite element mesh $\mathcal{T}_h(\theta_k)$ of the domain $\Omega(\theta_k)$ and an associated finite element approximation space $U_h(\theta_k) \subset U(\theta_k)$ defined by $U_h(\theta_k) = \{v_h(x) = \sum_{i=1}^{N(\theta_k)} \varphi_i(x, \theta_k)v_i, v_i \in \mathbb{R}\}$. The mesh, and then the basis functions $\varphi_i(x, \theta_k) \in U_h(\theta_k)$, clearly depend on the elementary event $\theta_k$.

The first drawback of these techniques is that they require numerous deterministic computations, each of which requiring the following classical steps: remeshing of the domain, assembling and solving a finite element problem of type (14). Another main drawback is that outcomes of the solution $u_h(\cdot, \theta_k) \in U_h(\theta_k)$ are defined on different meshes and that corresponding finite element problems (14) have a different algebraic structure. Information on non-computed outcomes is then unreachable. Moreover, an *a posteriori* probabilistic characterization of the solution is generally prohibitive since it would require storing the meshes and finite element solutions for all computed elementary events.

Of course, it is possible to focus on a mesh and geometry independent quantity of interest. Let us denote by $J(u_h(\cdot, \theta))$ this quantity of interest and by $J(u_h(\cdot, \theta_k))$ its computed outcomes. In the context of Monte-Carlo simulations, these outcomes allow evaluating statistics of $J$. In the case of projection methods (see section 3.3.2), outcomes correspond to stochastic integration
points which are used to compute the decomposition of J on a stochastic basis: $J = \sum_{\alpha \in \mathcal{J}} J_{\alpha} H_{\alpha}(\theta)$, with

$$J_{\alpha} = E(H_{\alpha}J(\mathbf{u}_h)) \approx \sum_{k} \omega_k H_{\alpha}(\theta_k) J(\mathbf{u}_h(\cdot, \theta_k)). \quad (25)$$

The problem is that these quantities of interest are not necessarily known a priori, before having an idea of the solution.

## 4.2 Finite element approximation on a random mesh

In order to avoid remeshings and to obtain an “explicit” description of the stochastic solution, an alternative technique consists in building a classical finite element approximation on a reference deterministic domain, and to introduce a random mapping between this domain and the random domain $\Omega(\theta)$. Such a strategy has been used in [13,14]. Let us briefly describe this methodology in our context to point out its advantages and drawbacks.

### 4.2.1 Description of the random geometry by a random mapping

We first introduce a reference deterministic domain $\Omega^0 \subset \mathbb{R}^d$ and suppose that there exists a continuous one to one mapping $T(\cdot, \theta)$ between $\Omega^0$ and the outcome $\Omega(\theta)$ of the random domain: $T(\cdot, \theta) : \mathbf{x}^0 \in \Omega^0 \mapsto \mathbf{x} = T(\mathbf{x}^0, \theta) \in \Omega(\theta)$. We denote by $T^*(\cdot, \theta)$ the inverse mapping. Reference domain $\Omega^0$ and mapping $T$ allow the description of the random domain $\Omega$. We suppose that the parts $\Gamma_i(\theta)$ of boundary $\partial\Omega(\theta)$ are the images by mapping $T(\cdot, \theta)$ of complementary parts $\Gamma_i^0$ of $\partial\Omega^0$. We denote by $T|_i$ the restriction of $T$ on $\Gamma_i^0$.

### 4.2.2 Construction of an approximation space

We introduce a fixed finite element mesh $\mathcal{T}_h^0$ of reference domain $\Omega^0$, leading to the definition of the following finite element approximation space: $\mathcal{U}_h^0 = \{\mathbf{v}^0(\mathbf{x}^0) = \sum_{i=1}^N \varphi_i^0(\mathbf{x}^0)v_i, v_i \in \mathbb{R}\}$, where the $\varphi_i^0 \in \mathcal{U}^0 = \{\mathbf{v} \in (H^1(\Omega^0))^d; \mathbf{v} = 0 \text{ on } \Gamma_1^0\}$ are the finite element basis functions associated with the mesh $\mathcal{T}_h^0$. Approximation space $\mathcal{U}_h^0$ being independent of the event, a tensor product approximation space can be introduced: $\mathcal{V}_{h,P}^0 \cong \mathbb{S}_P \otimes \mathcal{U}_h^0$. An approximation subspace $\mathcal{V}_{h,P}$ of $\mathcal{V}$ is then simply defined by a suitable mapping of functions in $\mathcal{V}_{h,P}^0$:

$$\begin{aligned} \mathcal{V}_{h,P} &= \{\mathbf{v}(\mathbf{x}, \theta) = \mathbf{v}^0(\mathbf{T}^*(\mathbf{x}, \theta), \theta), \mathbf{v}^0 \in \mathcal{V}_{h,P}^0\} \\ &= \{\mathbf{v}(\mathbf{x}, \theta) = \sum_{i=1}^N \varphi_i(\mathbf{x}, \theta)v_i(\theta), v_i(\theta) \in \mathbb{S}_P\}, \end{aligned} \quad (26)$$
where $\varphi_i(x, \theta) = \varphi_i^0(T^*(x, \theta))$. The image of $\mathcal{T}_h^0$ by mapping $T(\cdot, \theta)$ is a mesh of $\Omega(\theta)$, denoted by $\mathcal{T}_h(\theta)$, which can be considered as a random mesh.

### 4.2.3 Definition and construction of the approximation

At this point, different definitions of the approximation can still be used, as introduced in section 3.3. The approximation $u_{h,P} \in \mathcal{V}_{h,P}$ can be written $u_{h,P} = \varphi(x, \theta)u(\theta)$, with $u(\theta) = \sum_{\alpha \in \mathcal{J}_P} u_\alpha H_\alpha(\theta)$. The standard Galerkin approximation is still defined by problem (19). Although the approximation has not a tensor product structure, coefficients $u_\alpha \in \mathbb{R}^N$ are still defined by the system of equations (20), where components of random matrix $A$ and random vector $b$ can be written as integral on the reference domain. For every function $f(x, \theta)$, with $x \in \Omega(\theta)$, let $f^0(x^0, \theta) = f(T(x^0, \theta), \theta)$ be the corresponding function defined on $\Omega^0 \times \Theta$. Then, components of $A$ and $b$ write:

$$ (\mathbf{A}(\theta))_{ij} = a(\varphi_j, \varphi_i; \theta) = \int_{\Omega^0} \varepsilon(\varphi_i^0) : \mathbf{C}^0 : \varepsilon(\varphi_j^0) J^0 dx^0, \quad (27) $$

$$ (\mathbf{b}(\theta))_i = b(\varphi_i; \theta) = \int_{\Omega^0} \varphi_i^0 \cdot \mathbf{f}^0 J^0 dx^0 + \int_{\Gamma_2^0} \varphi_i^0 \cdot \mathbf{F}^0 J_2^0 ds^0, \quad (28) $$

where $J^0$ and $J_2^0$ are the jacobians of mappings $T$ and $T|_2$ respectively. In equation (27), the strain term can be written as follows: $\varepsilon(\varphi_j^0) = \frac{1}{2}(\frac{\partial \varphi_j^0}{\partial x} + \frac{\partial \varphi_j^0}{\partial x})$, with $\frac{\partial \varphi_j^0}{\partial x} = \frac{\partial \varphi_j^0}{\partial x^0} \cdot (\frac{\partial T}{\partial x})^{-1}$. The randomness on the geometry is then only contained in mappings $T$ and $T|_2$ and their jacobians. However, the integration of the left and right-hand sides of system (20) is not so trivial and necessitates a particular attention. An a priori explicit stochastic representation of the mappings can facilitate this integration.

Another definition of the approximation, based on a projection method, can also be used. It still consists in defining the random vector $u \in \mathbb{R}^N \otimes \mathcal{S}_P$ by equations (23) and (24). However, in this case, the approximation $u_{h,P}$ appears to be the projection of the semi-discretized solution with respect to the following inner product on $\mathcal{V}$:

$$ \begin{aligned} \ll \mathbf{u}, \mathbf{v} \gg &= \int_{\Theta} \int_{\Omega(\theta)} \mathbf{u} \cdot \mathbf{v} (J^0)^{-1} dx dP(\theta) \\ &= \int_{\Theta} \int_{\Omega^0} \mathbf{u}^0 \cdot \mathbf{v}^0 dx dP(\theta). \end{aligned} \quad (29) $$

### 4.2.4 Remarks

Compared to classical finite element techniques with remeshings, this method has the advantage to do not require any remeshing. It also leads to a complete description of the stochastic solution. However, some drawbacks can be pointed
out. First, the construction of mapping $T$ is not so trivial and necessitates the elaboration of ad-hoc numerical techniques [40,13]. Regularity requirements on the random mapping avoid dealing with general random domains, and in particular with topology changes. Moreover, since the approximation space has not a tensor product structure, a simple post-processing of random vector $u$ does not give meaningful information on the solution. For example, the expectation $E(u)$ leads to the expectation of $u^0(x^0, \theta)$ but not the expectation of the physical displacement field $u(x, \theta)$ at a given point $x$.

# 5 The eXtended Stochastic Finite Element Method (X-SFEM)

In this section, we present the basic principles of the eXtended Stochastic Finite Element Method (X-SFEM), which is the extension of the X-FEM method in the stochastic framework. In this article, we only deal with random shapes. This is a degenerated case of the X-FEM method [19,20] where only a classical finite element approximation can be used, without enrichment of the approximation space by the partition of unity method. In this section, the method is introduced from a theoretical point of view. Computational aspects related to the construction of the approximation will be detailed in the following section.

## 5.1 Reformulation of the problem on a deterministic domain

We introduce a deterministic spatial domain $B$ which contains all outcomes of the random domain, i.e. $\bigcup_{\theta \in \Theta} \Omega(\theta) \subset B$ (figure 3).

<figure><img src="image_4.png" /><figcaption>Fig. 3. Deterministic domain $B$ including all outcomes of the random domain $\Omega(\theta)$</figcaption></figure>

The problem is then reformulated by considering prolongation of functions of $\mathcal{V}$ on $B \times \Theta$. The aim of working in this deterministic domain is to facilitate the construction of approximation spaces, independently of the random geometry. We first introduce the function space

$$ \bar{\mathcal{U}}(\theta) = \{v \in (H^1(B))^d; v = 0 \text{ on } \Gamma_1(\theta)\}. \quad (30) $$
Then, the prolongation of the stochastic solution **u** will be searched in the following function space

$$
\begin{align}
\bar{\mathcal{V}} &= \{ \mathbf{v} : \theta \in \Theta \rightarrow \mathbf{v}(\cdot, \theta) \in \bar{\mathcal{U}}(\theta); \int_{\Theta} \|\mathbf{v}\|_{\bar{\mathcal{U}}(\theta)}^2 dP(\theta) < \infty \} \nonumber \\
&:= L^2(\Theta, dP; \bar{\mathcal{U}}) \tag{31}
\end{align}
$$

The stochastic problem (6) can now be equivalently reformulated as follows:
find **u** ∈ $\bar{\mathcal{V}}$ such that

$$
A(\mathbf{u}, \mathbf{v}) = L(\mathbf{v}) \quad \forall \mathbf{v} \in \bar{\mathcal{V}} \qquad (32)
$$

where bilinear form $A$ and linear form $L$ are defined as in the initial formulation (6). Bilinear form $A$ is clearly only semi-coercive on $\bar{\mathcal{V}}$, such that there exists infinitely many solutions of problem (32). However, we can show that the “physical part” of these solutions is unique and coincides with the solution of the initial problem (6) (see [15]). The physical part of a solution $\mathbf{u} \in \bar{\mathcal{V}}$ is the restriction $\mathbf{u}|_{\mathcal{P}}$ of $\mathbf{u}$ on the physical domain $\mathcal{P}$ of $B \times \Theta$, defined by

$$
\mathcal{P} = \{(\mathbf{x}, \theta) \in B \times \Theta; \mathbf{x} \in \Omega(\theta)\}. \quad (33)
$$

In other words, all solutions of problem (32) differs by functions whose support
is included in the complementary subset of P, which can be called the “non-
physical domain”. In the following, we will denote by I<sub>P</sub>(<b>x</b>, θ) the indicator
function of P, defined on B × Θ, which writes

$$
I_{\mathcal{P}}(\mathbf{x}, \theta) = \begin{cases} 1 & \text{if } (\mathbf{x}, \theta) \in \mathcal{P} \\ 0 & \text{otherwise} \end{cases} := I_{\Omega(\theta)}(\mathbf{x}), \qquad (34)
$$

where $I_{\Omega(\theta)}(\mathbf{x})$ is the classical indicator function of spatial domain $\Omega(\theta)$. Bi-linear form $A$ and linear form $L$ can be equivalently rewritten as follows:

$$
A(\mathbf{u}, \mathbf{v}) = \int_{\Theta} \int_{B} \varepsilon(\mathbf{v}) : \mathbf{C} : \varepsilon(\mathbf{u}) I_{\mathcal{P}} dx \, dP(\theta), \quad (35)
$$

$$
L(\mathbf{v}) = \int_{\Theta} \int_{B} \mathbf{v} \cdot \mathbf{f} \, I_{\mathcal{P}} \, dx \, dP(\theta) + \int_{\Theta} \int_{\Gamma_2(\theta)} \mathbf{v} \cdot \mathbf{F} \, ds \, dP(\theta). \quad (36)
$$

5.2 Representation of the geometry using random level-sets

5.2.1 Random level-sets

Here, we represent the random geometry in an implicit manner using the level-
set technique [18]. The level-set technique consists in representing an hyper-
surface $\Gamma$ in $\mathbb{R}^d$ by the iso-zero of a function $\phi$ called a level-set function. Then,
to represent a random hyper-surface $\Gamma(\theta)$, it is natural to introduce a random level-set function $\phi(\boldsymbol{x}, \theta)$ such that

$$ \Gamma(\theta) = \{\boldsymbol{x} \in \mathbb{R}^d; \phi(\boldsymbol{x}, \theta) = 0\}. \qquad (37) $$

The level-set is generally taken as the signed distance function to the hyper-surface: $\phi(\boldsymbol{x}, \theta) = \pm dist(\boldsymbol{x}, \Gamma(\theta))$. The distance function have some interesting properties. For example, on the hyper-surface, the gradient of $\phi$ with respect to $\boldsymbol{x}$ is the unitary normal to the hyper-surface. Let us note that the dirac measure on $\Gamma(\theta)$, characterized by a level-set $\phi$, can be written $\delta_{\Gamma(\theta)}(\boldsymbol{x}) = \delta(\phi(\boldsymbol{x}, \theta))$, where $\delta$ is the classical dirac delta function.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Random level-set $\phi(\boldsymbol{x}, \theta)$ whose iso-zero represents the boundary of random domain $\Omega(\theta)$</figcaption></figure>

Here, the hyper-surface to be represented can be the complete boundary $\partial\Omega(\theta)$ of domain $\Omega$ (figure 4) or a part of this boundary. By convention, we consider that $\phi(\cdot, \theta)$ takes negative values on random domain $\Omega(\theta)$ and positive values on the complementary domain in $B$. The random domain $\Omega(\theta)$ and its boundary can then be characterized by

$$ \Omega(\theta) = \{\boldsymbol{x} \in B; \phi(\boldsymbol{x}, \theta) < 0\}, \qquad (38) $$

$$ \partial\Omega(\theta) = \{\boldsymbol{x} \in B; \phi(\boldsymbol{x}, \theta) = 0\} \cup \{\boldsymbol{x} \in \partial B; \phi(\boldsymbol{x}, \theta) < 0\}. \qquad (39) $$

The indicator function of the physical domain, defined by (34), can then be simply written in terms of the level-set:

$$ I_{\mathcal{P}}(\boldsymbol{x}, \theta) = I_{\Omega(\theta)}(\boldsymbol{x}) = \mathcal{H}(-\phi(\boldsymbol{x}, \theta)) \qquad (40) $$

where $\mathcal{H}: \mathbb{R} \to \{0, 1\}$ is the heaviside function defined by $\mathcal{H}(y) = \begin{cases} 1 & \text{if } y > 0 \\ 0 & \text{if } y \le 0 \end{cases}$.

In order to define boundary conditions, we also need to describe random parts of the boundary $\partial\Omega$. This description is also greatly simplified by using the level-set technique. Let us see for example how to describe $\Gamma_2(\theta)$, where surface loads are applied, and how to define integrals on $\Gamma_2$. Two first cases can be encountered:
$$ \bullet \Gamma_2 = \{ \boldsymbol{x} \in B; \phi = 0 \} : $$

$$ \int_{\Gamma_2(\theta)} f \, ds = \int_{B} f \delta(\phi) \, dx \qquad (41) $$

$$ \bullet \Gamma_2 = \{ \boldsymbol{x} \in \partial_2 B; \phi < 0 \} : $$

$$ \int_{\Gamma_2(\theta)} f \, ds = \int_{\partial_2 B} f \mathcal{H}(-\phi) \, ds \qquad (42) $$

In (41), $\Gamma_2$ is the iso-zero of the level-set $\phi$. In (42), $\Gamma_2 = \partial_2 B \cap \Omega$, where $\partial_2 B$ is a deterministic part of the boundary of $B$. Of course, these cases do not span all the possibilities. For example, a specific treatment is required if $\Gamma_2$ is only included in the hyper-surfaces described in (41) and (42). A possibility consists in introducing an auxiliary level-set $\tilde{\phi}$ whose iso-zero characterizes the boundary of the hyper-surface. If we consider by convention that $\tilde{\phi}$ takes negative values on $\Gamma_2$, that leads to the following two other cases:

$$ \bullet \Gamma_2 = \{ \boldsymbol{x} \in B; \phi = 0; \tilde{\phi} < 0 \} : $$

$$ \int_{\Gamma_2(\theta)} f \, ds = \int_{B} f \mathcal{H}(-\tilde{\phi}) \delta(\phi) \, dx \qquad (43) $$

$$ \bullet \Gamma_2 = \{ \boldsymbol{x} \in \partial_2 B; \phi < 0; \tilde{\phi} < 0 \} : $$

$$ \int_{\Gamma_2(\theta)} f \, ds = \int_{\partial_2 B} f \mathcal{H}(-\phi) \mathcal{H}(-\tilde{\phi}) \, ds \qquad (44) $$

The case of equation (43) is illustrated on figure 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Level-set $\tilde{\phi}(\cdot, \theta)$ whose iso zero on $\partial\Omega(\theta)$ represents the boundary of hyper-surface $\Gamma_2(\theta)$</figcaption></figure>

### 5.2.2 On the construction of levels-sets

For some particular shapes (circles, hyper-planes, ellipses, polygons, etc), an explicit expression of level-sets can be used [20]. Then, considering the parameters of these level-sets as random variables (e.g. center coordinates and radius of the circle, etc), we easily define the associated random level-sets.
In the context of the level-set technique, a simple way to obtain more complex geometries consists in using classical boolean operations on basic domains. Let us suppose that $\Omega_1$ and $\Omega_2$ are two random domains, respectively characterized by random level-sets $\phi_1$ and $\phi_2$. Then, we want to define a random domain $\Omega$, characterized by a random level-set $\phi$, by doing boolean operations on $\Omega_1$ and $\Omega_2$. Below, we present classical boolean operations and the associated manipulations on level-sets:

* union: $\Omega = \Omega_1 \cup \Omega_2$,

$$ \phi(\mathbf{x}, \theta) = \min\{\phi_1(\mathbf{x}, \theta), \phi_2(\mathbf{x}, \theta)\} \quad (45) $$

* intersection: $\Omega = \Omega_1 \cap \Omega_2$,

$$ \phi(\mathbf{x}, \theta) = \max\{\phi_1(\mathbf{x}, \theta), \phi_2(\mathbf{x}, \theta)\} \quad (46) $$

* complement: $\Omega = \Omega_1^c$,

$$ \phi(\mathbf{x}, \theta) = -\phi_1(\mathbf{x}, \theta) \quad (47) $$

In this article, we will consider that level-sets are obtained by the above techniques and that the probabilistic characterization of these level-sets is given. The question of the probabilistic identification of random level-sets, as well as the representation of more general random shapes by level-sets, will be introduced in a subsequent paper.

### 5.2.3 Discretization of level-sets

In practice, a level-set is approximated at the space level by introducing a finite element mesh of domain $B$. This mesh can be different from the one which is used for the approximate solution. However, the use of the same finite element mesh greatly simplifies the computational aspects within the context of X-SFEM. We denote by $\{\varphi_i(\mathbf{x})\}_{i \in I}$ the set of finite element interpolation functions associated with the nodes $\{\mathbf{x}_i\}_{i \in I}$ of the mesh. A level-set, approximated at the space level, writes $\phi(\mathbf{x}, \theta) = \sum_{i \in I} \varphi_i(\mathbf{x})\phi_i(\theta)$, where the $\phi_i(\theta)$ are the nodal values of the random level-set, which are random variables. In this article, we will consider a linear interpolation for the level-set. The iso-zero of the approximate level-sets is then continuous and piecewise linear, thus approximating the geometry piecewise linearly.

The randomness on the geometry is completely contained in the random level-sets. Level-sets, after the stochastic modeling step (see section 2.4), can be expressed in terms of basic random variables $\xi$. In some cases, a simple and explicit representation of the level-set in terms of $\xi$ could be available. However, in practice, an approximation of the level-set at the stochastic level can
be convenient from a computational point of view. It consists in decomposing the level-set on a suitable stochastic basis (in practice the same as for the approximate solution): $\phi(\boldsymbol{x}, \boldsymbol{\xi}(\theta)) = \sum_{i \in I} \sum_{\alpha} \varphi_i(\boldsymbol{x}) H_{\alpha}(\boldsymbol{\xi}(\theta)) \phi_{i, \alpha}$. Giving the coefficients $\phi_{i, \alpha}$ then allow a complete characterization of the approximate random geometry and a quick evaluation of outcomes of the random level-set.

### 5.3 Discussion on boundary conditions

In the case where only Neumann boundary conditions are applied ($\Gamma_1 = \emptyset$), $\bar{\nabla} = L^2(\Theta, dP; (H^1(B))^d)$ has a tensor product structure, i.e. $\bar{\nabla} \cong L^2(\Theta, dP) \otimes (H^1(B))^d$. A classical tensor product approximation can then be easily constructed.

In the case where Dirichlet boundary conditions are applied on the random boundary's part $\Gamma_1$, $\bar{\nabla}$ has no more a tensor product structure. A reformulation of the problem must then be introduced if we want to work in a tensor product function space for the displacement field and then facilitates the construction of approximation spaces. A possible way to achieve this is to introduce a Lagrange multiplier $\boldsymbol{p}$ in order to impose the Dirichlet boundary conditions in a weak sense. The problem can be reformulated as follows: find $\boldsymbol{u} \in S \otimes (H^1(B))^d$ and $\boldsymbol{p} \in Q$ such that $\forall \boldsymbol{v} \in S \otimes (H^1(B))^d$ and $\forall \boldsymbol{q} \in Q$,

$$A(\boldsymbol{u}, \boldsymbol{v}) + \int_{\Theta} \int_{\Gamma_1(\theta)} \boldsymbol{p} \cdot \boldsymbol{v} \, ds \, dP(\theta) = L(\boldsymbol{v}), \quad (48)$$

$$\int_{\Theta} \int_{\Gamma_1(\theta)} \boldsymbol{q} \cdot \boldsymbol{u} \, ds \, dP(\theta) = 0. \quad (49)$$

A natural choice consists in introducing a boundary supported Lagrange multiplier, whose function space is defined by

$$\begin{aligned} Q = \{ \boldsymbol{p} : & \theta \mapsto \boldsymbol{p}(\cdot, \theta) \in H^{-1/2}(\Gamma_1(\theta)); \\ & \int_{\Theta} \|\boldsymbol{p}(\cdot, \theta)\|_{H^{-1/2}(\Gamma_1(\theta))}^2 dP(\theta) < \infty \}. \end{aligned} \quad (50)$$

Although the construction of an approximation space for the displacement field is then greatly simplified, a particular care must still be taken for the choice of the approximation space of $Q$ in order to satisfy the discrete inf-sup condition [41]. In the case where $\Gamma_1$ is deterministic, $Q \cong S \otimes H^{-1/2}(\Gamma_1)$ has a tensor product structure. The construction of approximation spaces can then be made independently at space and stochastic levels. It allows taking part of previous works in the deterministic context for a suitable construction of approximation subspaces of $(H^1(B))^d$ and $H^{-1/2}(\Gamma_1)$. For example, the methodology proposed in [42,43] can be used.

Of course, when $\Gamma_1$ is random, $Q$ has no more a tensor product structure. An alternative could consist in introducing a distributed Lagrange multiplier, with
$\mathcal{Q} \subset L^2(B) \otimes \mathcal{S}$, while keeping the same weak formulation. This allows to obtain a formulation with tensor product function spaces for both displacement field and Lagrange multiplier. Finally, to avoid the use of Lagrange multiplier, other alternatives could also be used to impose Dirichlet boundary conditions, such as penalty method [44] or Nitsche's method [45]. Those different strategies will be investigated in a subsequent article.

## 5.4 Definition of approximation spaces

From now on, we will consider that $\Gamma_1$ is deterministic and is contained in the boundary $\partial B$ of the deterministic domain (figure 6). The first condition makes that $\bar{\mathcal{V}}$ recovers a tensor product structure

$$ \bar{\mathcal{V}} = L^2(\Theta, dP; \bar{\mathcal{U}}) \cong L^2(\Theta, dP) \otimes \bar{\mathcal{U}}. \quad (51) $$

The second condition makes classical the construction of a finite element subspace in $\bar{\mathcal{U}}$, which satisfies the Dirichlet boundary conditions.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Case of a deterministic boundary part $\Gamma_1 \subset \partial B$</figcaption></figure>

Let $\mathcal{T}_h$ be a finite element mesh of domain $B$ and $\{\varphi_i\}_{i=1}^N$ be a finite element approximation basis in $\bar{\mathcal{U}}$. The approximation space at the space level can then be written:

$$ \bar{\mathcal{U}}_h = \{ \mathbf{v}(\mathbf{x}) = \sum_{i=1}^{N} \varphi_i(\mathbf{x}) v_i = \varphi(\mathbf{x}) \mathbf{v}, \mathbf{v} \in \mathbb{R}^N \}. \quad (52) $$

A semi-discretized approximation space $\bar{\mathcal{V}}_h \subset \bar{\mathcal{V}}$ is defined as follows:

$$ \bar{\mathcal{V}}_h = \mathcal{S} \otimes \bar{\mathcal{U}}_h = \{ \mathbf{v}(\mathbf{x}, \theta) = \varphi(\mathbf{x}) \mathbf{v}(\theta), \mathbf{v} \in \mathcal{S} \otimes \mathbb{R}^N \}. \quad (53) $$

At the stochastic level, we introduce the approximation space $\mathcal{S}_P \subset \mathcal{S}$, defined in (17). Finally, as in equation (18), we introduce the following approximation space in $\bar{\mathcal{V}}$:

$$ \bar{\mathcal{V}}_{h,P} = \mathcal{S}_P \otimes \bar{\mathcal{U}}_h = \{ \mathbf{v}(\mathbf{x}, \theta) = \varphi(\mathbf{x}) \mathbf{v}(\theta), \mathbf{v} \in \mathcal{S}_P \otimes \mathbb{R}^N \}. \quad (54) $$
## 5.5 Galerkin approximation

The Galerkin approximation $\mathbf{u}_{h,P} \in \bar{\mathcal{V}}_{h,P}$ is defined by

$$A(\mathbf{u}_{h,P}, \mathbf{v}_{h,P}) = L(\mathbf{v}_{h,P}) \quad \forall \mathbf{v}_{h,P} \in \bar{\mathcal{V}}_{h,P}, \qquad (55)$$

which is equivalent to the following system of equations:

$$\sum_{\alpha \in \mathcal{J}_P} E(H_{\beta} A H_{\alpha}) \mathbf{u}_{\alpha} = E(b H_{\beta}) \quad \forall \beta \in \mathcal{J}_P. \qquad (56)$$

Computational aspects related to the resolution of problem (56) will be detailed in the following section. Below, we show that the Galerkin approximation can be said optimal in a certain sense with respect to the exact solution. Since A is symmetric, continuous and coercive on $\mathcal{V}$, it defines an inner product on $\mathcal{V}$, defined $\forall \mathbf{u}, \mathbf{v} \in \mathcal{V}$ by

$$\ll \mathbf{u}, \mathbf{v} \gg_A = A(\mathbf{u}, \mathbf{v}). \qquad (57)$$

We denote by $\| \cdot \|_A$ the associated norm. By prolongation, $\ll$, $\gg_A$ (resp. $\| \cdot \|_A$) defines a semi-inner product (resp. semi-norm) on $\bar{\mathcal{V}}$. Let us denote by $\mathbf{u} \in \bar{\mathcal{V}}$ a solution of problem (32). Then, the Galerkin approximation $\mathbf{u}_{h,P}$ is the optimal approximation of $\mathbf{u}$ in $\bar{\mathcal{V}}_{h,P}$ with respect to the semi-norm $\| \cdot \|_A$:

$$\| \mathbf{u}_{h,P} - \mathbf{u} \|_A = \min_{\mathbf{v}_{h,P} \in \bar{\mathcal{V}}_{h,P}} \| \mathbf{v}_{h,P} - \mathbf{u} \|_A. \qquad (58)$$

Equation (58) can also be interpreted as follows: $\mathbf{u}_{h,P}$ is the projection of $\mathbf{u}$ on $\bar{\mathcal{V}}_{h,P}$ with respect to the inner-product $\ll$, $\gg_A$. The semi-inner product giving no weight to the non-physical part of the solution, it can also be interpreted as follows: the physical part of the approximate solution $\mathbf{u}_{h,P|\mathcal{P}}$ is the projection of the exact physical solution $\mathbf{u}_{|\mathcal{P}}$ with respect to the inner product $\ll$, $\gg_A$ on $\mathcal{V}$. We will see in examples that the Galerkin approximation leads to a very good approximation of the exact physical solution.

**Remark 2 (Non-uniqueness of the approximation)** System (56) may be only semi-definite, which means that the Galerkin approximation may not be unique. It is the case if a basis function $\varphi_i(\mathbf{x}) H_{\alpha}(\theta)$ is in the left kernel of bilinear form A, which is only semi-coercive on $\bar{\mathcal{V}}$. This occurs if the measure of the intersection between the support of the basis function and the physical domain $\mathcal{P}$ is zero, i.e.

$$\mu(\mathrm{supp}(\varphi_i(\mathbf{x}) H_{\alpha}(\theta)) \cap \mathcal{P}) = 0, \qquad (59)$$

where $\mu = \lambda \otimes P$ is the product measure on $B \times \Theta$, with $\lambda$ the Lebesgue measure on B. The support of the basis function is defined by

$$\mathrm{supp}(\varphi_i(\mathbf{x}) H_{\alpha}(\theta)) = \mathrm{closure}(\{(\mathbf{x}, \theta) \in B \times \Theta; \varphi_i(\mathbf{x}) H_{\alpha}(\theta) \neq 0\}). \qquad (60)$$
The associated degree of freedom is undetermined but since it only contributes to the non-physical part of the solution, it can be arbitrarily chosen. In practise, they can be easily detected in the resolution procedure since they correspond to zero diagonal terms of the system of equation (56), i.e. $E((\mathbf{A})_{ii}H_{\alpha}H_{\alpha}) = 0$. Of course, one can also use any linear solver allowing the obtention of a particular solution of a singular linear system of equations.

In the particular case where $\text{supp}(H_{\alpha}) = \Theta$, which is the case when we use classical polynomial chaos basis, condition (59) is equivalent to $P(\{\lambda(\text{supp}(\varphi_i) \cap \Omega) = 0\}) = 1$ and occurs if the set of elements in $\mathcal{T}_h$ composing $\text{supp}(\varphi_i)$ is almost surely outside the domain $\Omega$. Elements surely outside the domain being automatically detected by the resolution procedure proposed in the following section, the detection of undetermined degrees of freedom is straightforward. In the contrary, if all elements in $\mathcal{T}_h$ have a non-zero probability to intersect the domain $\Omega$, the Galerkin approximation is unique.

## 5.6 About projection methods

### 5.6.1 A $L^2$ projection method based on X-FEM ($P$-$X$-FEM)

Here, we present another definition of the approximation which can be obtained by simply applying a classical $L^2$ projection method (see section 3.3.2), where the deterministic code is based on the deterministic X-FEM method. Let $\mathbf{u}_h \in \bar{\mathcal{V}}_h$ be a solution of the semi-discretized problem:

$$A(\mathbf{u}_h, \mathbf{v}_h) = L(\mathbf{v}_h) \quad \forall \mathbf{v}_h \in \bar{\mathcal{V}}_h. \qquad (61)$$

Denoting by $\mathbf{u}_h = \varphi(\mathbf{x})\mathbf{u}(\theta)$, the random vector $\mathbf{u}(\theta)$ verifies P-almost surely

$$\mathbf{A}(\theta)\mathbf{u}(\theta) = \mathbf{b}(\theta). \qquad (62)$$

Equation (62) corresponds to the classical system of equations which could be obtained by the X-FEM method, when considering the deterministic problem associated with elementary event $\theta$. $\mathbf{A}(\theta)$ may be only semi-definite positive. A natural definition of a pseudo-inverse $\mathbf{A}(\theta)^{+}$ of $\mathbf{A}(\theta)$ consists in setting to zero the undetermined degrees of freedom. They correspond to basis functions $\varphi_i$ whose support is such that $\lambda(\text{supp}(\varphi_i) \cap \Omega(\theta)) = 0$, where $\lambda$ is the classical Lebesgue measure, or equivalently to basis functions $\varphi_i$ such that all finite elements composing $\text{supp}(\varphi_i)$ are outside the domain $\Omega(\theta)$.

A classical projection method, as introduced in section 3.3.2, would then consist in projecting $\mathbf{u}$ on $S_P \otimes \mathbb{R}^N$ with respect to the natural inner-product in $L^2(\Theta, dP; \mathbb{R}^N)$. The coefficients $\mathbf{u}_{\alpha}$, $\alpha \in \mathcal{J}_P$, of the approximate solution
$\mathbf{u}_{h,P} \in \bar{\mathcal{V}}_{h,P}$ are then approximated by a suitable stochastic integration

$$ \mathbf{u}_{\alpha} = E(H_{\alpha}\mathbf{u}) \approx \sum_{k} \omega_{k} H_{\alpha}(\theta_{k}) \mathbf{A}^{+}(\theta_{k})\mathbf{b}(\theta_{k}). \quad (63) $$

At the continuous level, this is equivalent to define $\mathbf{u}_{h,P}$ as the projection of $\mathbf{u}_h \in \bar{\mathcal{V}}_h$ on $\bar{\mathcal{V}}_{h,P}$ in the following sense:

$$ \|\mathbf{u}_{h,P} - \mathbf{u}_h\|_{L^2(B \times \Theta, d\mu)} = \min_{\mathbf{v}_{h,P} \in \bar{\mathcal{V}}_{h,P}} \|\mathbf{v}_{h,P} - \mathbf{u}_h\|_{L^2(B \times \Theta, d\mu)}, \quad (64) $$

where $\mu = \lambda \otimes P$ is the product measure on $B \times \Theta$, and

$$ \|\mathbf{v}\|_{L^2(B \times \Theta, d\mu)}^2 = \int_{\Theta} \int_{B} \mathbf{v}(\mathbf{x}, \theta) \cdot \mathbf{v}(\mathbf{x}, \theta) d\mathbf{x} d\mathbf{P}(\theta). \quad (65) $$

As we will see in examples, this strategy leads to a relatively bad approximation, compared to the Galerkin approximation. In fact, significant errors are made on finite elements $K \in \mathcal{T}_h$ which have a non-zero probability to be outside the domain but also a non-zero probability to intersect the domain, i.e. $0 < P(\{\lambda(K \cap \Omega(\theta)) = 0\}) < 1$.

### 5.6.2 What should be a good projection method based on X-FEM?

The problem with the above classical projection method is that it gives a weight to the non-physical part of the solution, obtained by an artificial prolongation to zero of the solution. In fact, a good projection technique should be based on a physical norm of the solution, as for the Galerkin approach. For example, we could define the approximate solution by

$$ \|\mathbf{u}_{h,P} - \mathbf{u}_h\|_{L^2(\mathcal{T}, d\mu)} = \min_{\mathbf{v}_{h,P} \in \bar{\mathcal{V}}_{h,P}} \|\mathbf{v}_{h,P} - \mathbf{u}_h\|_{L^2(\mathcal{T}, d\mu)}, \quad (66) $$

where $\mathbf{u}_h(\mathbf{x}, \theta) = \varphi(\mathbf{x})\mathbf{u}(\theta)$ is the semi-discretized solution and where the norm $\| \cdot \|_{L^2(\mathcal{T}, d\mu)}$ is a classical $L^2$-norm on the restriction of functions on the physical domain $\mathcal{T}$. That leads to solve the following stochastic problem:

$$ \int_{\Theta} \int_{\Omega(\theta)} \mathbf{v}_{h,P} \cdot \mathbf{u}_{h,P} d\mathbf{x} d\mathbf{P}(\theta) = \int_{\Theta} \int_{\Omega(\theta)} \mathbf{v}_{h,P} \cdot \mathbf{u}_h d\mathbf{x} d\mathbf{P}(\theta). \quad (67) $$

Problem (67) is equivalent to the following system of equations:

$$ \sum_{\alpha \in \mathcal{J}_P} E(H_{\beta} \mathbf{M} H_{\alpha}) \mathbf{u}_{\alpha} = E(H_{\beta} \mathbf{M} \mathbf{u}) \quad \forall \beta \in \mathcal{J}_P, \quad (68) $$
where $\mathbf{M}: \Theta \mapsto \mathbb{R}^{N \times N}$ is the random geometrical mass matrix, whose coefficients are defined by

$$ (\mathbf{M}(\theta))_{ij} = \int_{\Omega(\theta)} \varphi_i(\mathbf{x}) \cdot \varphi_j(\mathbf{x}) \, dx \\ = \int_B \varphi_i(\mathbf{x}) \cdot \varphi_j(\mathbf{x}) I_{\mathcal{P}}(\mathbf{x}, \theta) \, dx. \qquad (69) $$

A good projection method, based on a projection of the physical part of the solution, then leads to the resolution of a problem similar to (56), which is obtained by the Galerkin method. This “good” projection method based on X-FEM would lead to the same computational aspects as X-SFEM.

# 6 Computational aspects of X-SFEM

In this section, we focus on the computational aspects of the X-SFEM method, i.e. the construction and resolution of the discretized problem (55), leading to the Galerkin approximation. We will first present the general resolution procedure, which is quite classical in the context of stochastic finite element methods. Then, we will come back on numerical aspects which are specific to the X-SFEM method. From now on, we consider that the probabilistic content of the problem is represented by a finite set of random variables $\xi$. We denote by $(\Theta, \mathcal{B}, P_{\xi})$ the associate finite-dimensional probability space. All random quantities will then be expressed in terms of $\xi$ instead of $\theta$ (see section 2.4).

## 6.1 The general resolution procedure

The coefficients of the Galerkin approximation $\mathbf{u}_{h,P} \in \bar{\mathcal{V}}_{h,P}$ are solution of the following system of equations:

$$ \sum_{\alpha \in \mathcal{J}_P} E(H_{\beta} \mathbf{A} H_{\alpha}) \mathbf{u}_{\alpha} = E(\mathbf{b} H_{\beta}) \quad \forall \beta \in \mathcal{J}_P, \qquad (70) $$

where $\mathbf{A}$ and $\mathbf{b}$ are respectively the random finite element matrix and vector. A first step generally consists in decomposing these quantities on the orthonormal basis $\{H_{\alpha}\}_{\alpha \in \mathcal{J}}$ of $\mathcal{S} = L^2(\Theta, dP_{\xi})$:

$$ \mathbf{A} = \sum_{\gamma \in \mathcal{J}_A} \mathbf{A}_{\gamma} H_{\gamma} \quad \text{with} \quad \mathbf{A}_{\gamma} = E(\mathbf{A} H_{\gamma}), \qquad (71) $$

$$ \mathbf{b} = \sum_{\gamma \in \mathcal{J}_b} \mathbf{b}_{\gamma} H_{\gamma} \quad \text{with} \quad \mathbf{b}_{\gamma} = E(\mathbf{b} H_{\gamma}). \qquad (72) $$

From orthonormality of the basis functions, the coefficients of the decomposition of $\mathbf{b}$ lead directly to the right-hand side terms of system (70). Then,
these coefficients must be computed only for indices $J_b = J_P$. The left hand side terms of system (70) can be written

$$E(\mathbf{A}H_{\alpha}H_{\beta}) = \sum_{\gamma \in J_A} \mathbf{A}_{\gamma}E(H_{\gamma}H_{\alpha}H_{\beta}). \qquad (73)$$

When using polynomial (eventually piecewise polynomial) basis functions {$H_{\alpha}$}, if $J_P$ corresponds to polynomial up to degree $p$, the coefficients $\mathbf{A}_{\gamma}$ must only be computed for the set of indices $J_A$ corresponding to polynomial functions up to degree $2p$. Indeed, $E(H_{\gamma}H_{\beta}H_{\alpha}) = 0$ for $\gamma \in (J \setminus J_A)$. Computing $E(H_{\gamma}H_{\alpha}H_{\beta})$ is classical within the context of Galerkin stochastic finite element methods. These quantities only depend on the chosen stochastic basis and are generally pre-computed. In the case where we use classical polynomial chaos or generalized polynomial chaos, an analytical expression of these quantities is sometimes available.

Finally, system (70), which is a huge system of size $P \times N$, is classically solved by iterative solvers as preconditioned conjugate gradient [46,6,7]. These solvers take part of the sparsity of the system, coming both from the sparse structure of random matrix $\mathbf{A}$ and from orthogonality properties of basis functions in $S_P$. Some alternative resolution techniques can be used at this step for saving computational times and memory requirements [3,8].

## 6.2 Computing the stochastic decomposition of random matrix and vector

In the general resolution procedure presented above, the only step which requires further details concerns the stochastic decomposition of random matrix and vector. To simplify the presentation, we consider that material properties (represented by tensor $\mathbf{C}$) and loadings ($\mathbf{F}$ and $\mathbf{f}$) are deterministic. Also to simplify the presentation, we consider that surface loads are applied on a random part of the boundary $\Gamma_2$ which can be characterized by the level-set $\phi$ as follows: $\Gamma_2(\xi) = \{\mathbf{x} \in B; \phi(\mathbf{x}, \xi) = 0\} \cup \{\mathbf{x} \in \partial_2 B; \phi(\mathbf{x}, \xi) < 0\}$, where $\partial_2 B$ is a part of the boundary $\partial B$. This is a combination of the two cases introduced in (41) and (42).

### 6.2.1 Decomposition of elementary matrices and vectors

Finite element matrix and vector are obtained by assembling element contributions:

$$\mathbf{A}(\xi) = \sum_{K \in \mathcal{J}_h} \mathbf{A}_K(\xi), \quad \mathbf{b}(\xi) = \sum_{K \in \mathcal{J}_h} \mathbf{b}_K(\xi), \qquad (74)$$
where components of elementary matrices and vectors are defined by

$$ (\mathbf{A}_K)_{ij} = \int_K \varepsilon(\varphi_j) : \mathbf{C} : \varepsilon(\varphi_i) \mathcal{H}(-\phi) \, dx, \qquad (75) $$

$$ (\mathbf{b}_K)_i = \int_K \varphi_i \cdot \mathbf{f} \mathcal{H}(-\phi) \, dx + \int_{\partial K \cap \partial_2 B} \varphi_i \cdot \mathbf{F} \mathcal{H}(-\phi) \, ds + \int_K \varphi_i \cdot \mathbf{F} \delta(\phi) \, dx. \qquad (76) $$

The stochastic decompositions (71) and (72) can be performed element by element. Then, we need to compute coefficients $\{\mathbf{A}_{K,\gamma}\}_{\gamma \in \mathcal{J}_A}$ and $\{\mathbf{b}_{K,\gamma}\}_{\gamma \in \mathcal{J}_b}$, defined by

$$ \mathbf{A}_{K,\gamma} = E(\mathbf{A}_K(\boldsymbol{\xi}) H_\gamma(\boldsymbol{\xi})) = \int_{\Theta} \mathbf{A}_K(\mathbf{y}) H_\gamma(\mathbf{y}) \, dP_{\boldsymbol{\xi}}(\mathbf{y}), \qquad (77) $$

$$ \mathbf{b}_{K,\gamma} = E(\mathbf{b}_K(\boldsymbol{\xi}) H_\gamma(\boldsymbol{\xi})) = \int_{\Theta} \mathbf{b}_K(\mathbf{y}) H_\gamma(\mathbf{y}) \, dP_{\boldsymbol{\xi}}(\mathbf{y}). \qquad (78) $$

Coefficients of the decomposition can be obtained by a suitable integration on $\Theta$ with respect to probability measure $dP_{\boldsymbol{\xi}}$. Let us denote by $Q_K = \{(\boldsymbol{\xi}, \omega)\}$ a quadature rule, where $\boldsymbol{\xi}$ and $\omega$ denote respectively integration points and integration weights, such that

$$ E(f(\boldsymbol{\xi})) = \int_{\Theta} f(\mathbf{y}) dP_{\boldsymbol{\xi}}(\mathbf{y}) \approx \sum_{(\boldsymbol{\xi}, \omega) \in Q_K} \omega f(\boldsymbol{\xi}). \qquad (79) $$

Computing coefficients of the decomposition then requires computing element quantities $\mathbf{A}_K(\boldsymbol{\xi})$ and $\mathbf{b}_K(\boldsymbol{\xi})$ for some given elementary events $\boldsymbol{\xi}$, corresponding to given outcomes $\phi(\cdot, \boldsymbol{\xi})$ of the random level-set. These computations are classical within the context of the deterministic X-FEM method [20]. If the element is cut by the iso-zero of the level set $\phi(\cdot, \boldsymbol{\xi})$, a suitable subdivision of the element and of its edges, illustrated on figure 7, allows defining a spatial numerical integration on $K \cap \Omega(\boldsymbol{\xi})$, $\partial K \cap \Gamma_2(\boldsymbol{\xi})$ and $K \cap \Gamma_2(\boldsymbol{\xi})$.

The question is then: how to define a suitable quadrature rule $Q_K$ on $\Theta$? Of course, it depends on the regularity of $\mathbf{A}_K(\boldsymbol{\xi})$ and $\mathbf{b}_K(\boldsymbol{\xi})$ with respect to $\boldsymbol{\xi}$. For example, if components of $\mathbf{A}_K(\boldsymbol{\xi})$ were polynomial of degree up to $j$ with respect to $\boldsymbol{\xi}$, and if we were using polynomial basis functions $H_\alpha$ up to degree $p$, the integration of $(\mathbf{A}_K H_\alpha)$ on $\Theta$ would require the integration of polynomial functions of degree $p + j$. We could then construct a classical Gaussian quadrature rule associated with measure $dP_{\boldsymbol{\xi}}$, leading to an exact computation of coefficients $\mathbf{A}_{K,\gamma}$.
<figure><img src="image_8.png" /><figcaption>Fig. 7. For a given outcome $\phi(\cdot, \xi)$ of the level-set, subdivision of a 3-nodes triangle $K$ and of its edges for the numerical integration on $K \cap \Omega(\xi)$, $\partial K \cap \Gamma_2(\xi)$ and $K \cap \Gamma_2(\xi)$.</figcaption></figure>

### 6.2.2 Stochastic regularity of element matrices and vectors

Let us take a simple example to analyze the regularity of element matrices and vectors with respect to random variables $\xi$. We consider that the spatial mesh $\mathcal{T}_h$ is composed by 3-nodes triangles. We also consider that the material is homogeneous and deterministic. Basis functions $\varphi_i$ being linear on an element $K \in \mathcal{T}_h$, element matrix $\mathbf{A}_K(\xi)$ writes

$$ \mathbf{A}_K(\xi) = \mathbf{A}_K^0 \int_K \mathcal{H}(-\phi(\mathbf{x}, \xi)) \, dx = \mathbf{A}_K^0 \lambda(K \cap \Omega(\xi)), \quad (80) $$

where $(\mathbf{A}_K^0)_{ij} = \varepsilon(\varphi_i) : C : \varepsilon(\varphi_j)$ is a constant matrix. In this case, the stochastic regularity of $\mathbf{A}_K$ only depends on the regularity of the function $\lambda(K \cap \Omega(\xi))$, which is the measure of the intersection between $K$ and the physical domain. Let us consider that the iso-zero of the level-set is an hyper-plane with a fixed unitary normal $\mathbf{a}$ and a position which depends on a unique random variable $\xi$. The level-set can be written $\phi(\mathbf{x}, \xi) = \mathbf{x} \cdot \mathbf{a} - \xi$. Figure 8 illustrates the dependence on $\xi$ of the function $\lambda(K \cap \Omega(\xi))$. The stochastic domain, included in $\mathbb{R}$, can be split into four intervals. On $]-\infty, \xi_1[$, $\lambda(K \cap \Omega(\xi))$ is equal to zero, on $]\xi_3, \infty[$, it is constant and equal to $\lambda(K)$. On the two other intervals $]\xi_1, \xi_2[$ and $]\xi_2, \xi_3[$, $\lambda(K \cap \Omega(\xi))$ is polynomial of degree 2 with respect to $\xi$. These last two intervals correspond to the case where the iso-zero of the level-set is cutting two edges of the element, each interval being associated with a different pair of edges.

In the same way, if surface loads $\mathbf{F}$ and body loads $\mathbf{f}$ are polynomial functions of $\mathbf{x}$, we can easily show that the surface loads and body loads contributions in $\mathbf{b}_K$ are piecewise polynomial with respect to $\xi$. Then, it is clear that such functions of $\xi$ can't be integrated by a classical quadrature constructed on the whole stochastic domain. The idea is then to build a suitable partition of the stochastic domain, on which finite element quantities are piecewise regular, in order to define a suitable quadrature rule.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Function $\lambda(K \cap \Omega(\xi))$ as a function of $\xi$. Level-set $\phi(x, \xi) = x \cdot a - \xi$</figcaption></figure>

In fact, the above partition of $\Theta$ can be completely characterized by the values
of the random level-set. Indeed, we can observe that each stochastic subdomain
corresponds to a constant state of the element $K$, defined by the sign of the
level-set at the nodes of the element. For our example, the four states of the
element, corresponding to each subdomain of $\Theta$, are illustrated on figure 9.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Sign of the level-set $\phi(\cdot, \xi)$ at nodes of element $K$ for each stochastic subdomain.</figcaption></figure>

### 6.2.3 Partitioning the stochastic domain for integration

Following the above observations, we here introduce a general procedure to
create a partition of the stochastic domain which is adapted to a finite element
$K$. This will allow us to define a suitable quadrature rule on $\Theta$ for integrating
element matrix and vector. This procedure only requires the use of the set
of values of the level-set at the nodes $I_K$ of the finite element, denoted by $(\phi_i(\xi))_{i \in I_K}$.

We suppose that random variables $\xi = (\xi_1, \dots, \xi_m)$ are statistically independent and that each $\xi_j$ is a uniform random variable with values in $[a_j, b_j]$. This can be obtained by a suitable iso-probabilistic mapping of the basic random variables. Then, the stochastic domain is the hyper-rectangle $\Theta = \prod_{j=1}^m [a_j, b_j] \subset \mathbb{R}^m$.

For a given elementary event $\xi$, we define the state of an element by $S_K(\xi) = (\mathrm{sign}(\phi_i(\xi)))_{i \in I_K}$, where $\mathrm{sign}(\phi_i(\xi)) \in \{-1, 0, 1\}$. The aim is to create a partition $C_K$ of $\Theta$ such that on each subdomain $C \in C_K$, $S_K(\xi)$ takes the same value for all $\xi \in C$. The proposed procedure consists in building a $m$-dimensional cartesian mesh $C_K$ by a “tree type” recursive procedure. An element $C \in C_K$, called a cell, is an hyper-rectangle $C = \prod_{j=1}^m [\xi_j^{(1)}, \xi_j^{(2)}]$. It has $2^m$ vertices, denoted by $\mathrm{vertices}(C)$. We denote by $\mathrm{split}(C)$ the boolean function indicating whether a cell should be split or not. The boolean value $\mathrm{split}(C)$ is false if $S_K(\xi) = S_K(\xi')$ for all $\xi, \xi' \in \mathrm{vertices}(C)$, and true otherwise. We denote by $\mathrm{order}(C, j)$ the order of the cell along dimension $j$, defined as the number of splittings along dimension $j$ which has led to $C$, starting from $\Theta$. Then, we define the order of a cell by $\mathrm{order}(C) = \max_j(\mathrm{order}(C, j))$ and the order of a partition by $\mathrm{order}(C_K) = \max_{C \in C_K}(\mathrm{order}(C))$. Let us denote by $k$ the maximum order of the partition, which will be directly related to the precision of the associated quadrature rule.

When splitting a cell $C$, there are different ways of defining the children cells. The basic way consists in performing an isotropic splitting, leading to $2^m$ children. We denote by $\mathrm{children}(C; k)$ the set of children cells of $C$ with order less than $k$, which means that $\mathrm{children}(C; k)$ is empty if $\mathrm{order}(C) = k$. However, in order to reduce the number of cells in the partition and then to reduce the size of the quadrature rule on $\Theta$, it can be interesting to use an anisotropic splitting criterium, defined as follows. A cell $C$ has $2m$ faces, grouped by pairs $(\mathcal{F}_j, \mathcal{F}'_j)$ of opposite faces orthogonal to dimension $j$. A parent cell will then be split along dimension $j$ only if $\mathrm{order}(C, j) < k$ and if $\mathrm{state}(\mathcal{F}_j) \neq \mathrm{state}(\mathcal{F}'_j)$, where the state of a face $\mathcal{F}$ is defined by $\mathrm{state}(\mathcal{F}) = \{S_K(\xi)\}_{\xi \in \mathrm{vertices}(\mathcal{F})}$. Vertices of two opposite faces $\mathcal{F}$ and $\mathcal{F}'$ are ordered in the same way. Then, they have the same state if each pair $(\xi, \xi')$ of opposite vertices have the same state, i.e. $S_K(\xi) = S_K(\xi')$. We denote by $\mathrm{children}(C; k)$ the obtained set of children of order less than $k$. Let us denote that for a given cell, the anisotropic criterium can lead to different non-empty sets of children for different $k$ values.

The meshing algorithm then only consists in applying function `recursive_split`, described in algorithm 1, to the initial domain $\Theta$ or to cells of a pre-defined partition of $\Theta$.
**Algorithm 1** Recursive construction of a m-dimensional cartesian mesh $C_K$ of order k

function recursive_split(C)

1: **if** split(C) and children(C; k) not empty **then**

2: **for all** C' ∈ children(C; k) **do**

3: recursive_split(C')

4: **end for**

5: **else**

6: add C to mesh $C_K$

7: **end if**

Figures 10 and 11 illustrate the isotropic and anisotropic splitting respectively for a 2D and 3D cells.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison between isotropic and anisotropic split of a 2D cell. □ and △ denote two different states of the vertices.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison between isotropic and anisotropic split of a 3D cell. □ and △ denote two different states of the vertices.</figcaption></figure>

Finally, a quadrature rule $Q_K$ on $\Theta$ can be defined by using the resulting cartesian mesh $C_K$ of $\Theta$. For that, we proceed as follows. For each cell $C \in C_K$, let $f_C$ be the linear mapping which transforms the reference hyper-cube $\tilde{C} = [0, 1]^m$ into $C = \prod_{j=1}^m [\xi_j^{(1)}, \xi_j^{(2)}]$ and let $\tilde{Q}$ be a classical Gauss-Legendre quadrature rule on $\tilde{C}$. Let us denote by $(\tilde{\xi}, \tilde{\omega}) \in \tilde{Q}$ the corresponding integration points and weights. A quadrature rule $Q_{C,K}$ on $C$ is then defined as follows:

$$Q_{C,K} = \{(\xi, \omega) = (f_C(\xi), \tilde{\omega}P_\xi(C)), (\tilde{\xi}, \tilde{\omega}) \in \tilde{Q}\}, \text{ where } P_\xi(C) = \frac{\prod_{j=1}^m (\xi_j^{(2)} - \xi_j^{(1)})}{\prod_{j=1}^m (b_j - a_j)}$$
is the probability of $C$. Finally, $Q_K$ is defined by $Q_K = \cup_{C \in C_K} Q_{C,K}$.

### 6.2.4 Summary

The stochastic decomposition of element matrices $A_K$ and element vectors $b_K$ then requires the following steps:
• Creation of the partition $C_K$ of $\Theta$ with given maximum order $k$

• Definition of a quadrature rule $Q_K$ associated with $C_K$

• Computation of $\mathbf{A}_K(\xi)$ and $\mathbf{b}_K(\xi)$ for each integration point $\xi$ in $Q_K$ (spatial integration defined in equations (75) and (76))

• Computation of coefficients $\mathbf{A}_{K,\gamma}$ and $\mathbf{b}_{K,\gamma}$ (defined by equations (77) and (78)).

All these steps, which require the main computational effort of X-SFEM, must be performed for each finite element $K \in \mathcal{T}_h$. However, these calculations are independent of one another so that the stochastic integration can be completely parallelized at the finite element level.

**Remark 3** The quadrature rule on each subdomain of stochastic partitions will be chosen in accordance with the expected theoretical regularity of $\mathbf{A}_K$ and $\mathbf{b}_K$, in order to well integrate $\mathbf{A}_K(\xi)H_\alpha(\xi)$ and $\mathbf{b}_K(\xi)H_\alpha(\xi)$. This regularity should be estimated a priori. The regularity of $\mathbf{A}_K$ and $\mathbf{b}_K$ depends on the type of finite elements and of the regularity of the level-set with respect to $\xi$. For example, for the case of a 3-nodes triangle illustrated in section 6.2.2, it was found that $\mathbf{A}_K(\xi)$ was piecewise polynomial of degree 2. In this case, the order of the quadrature can then be chosen in order to exactly integrate polynomials of degree $(2+j)$, where $j$ is the maximum degree of polynomials $H_\alpha(\xi)$.

# 7 Numerical examples

In this section, three numerical examples will illustrate the efficiency of the X-SFEM method. In example 1, we consider a random plate in tension where a probabilistic modeling with two independent random variables is chosen. As this problem has an analytical solution, true error indicators will be introduced in order to estimate the precision of the X-SFEM method. In this first problem, the approximation space contains the exact solution, so that X-SFEM leads to the exact numerical solution. This will allow us to focus on the precision of the integration at the stochastic level. Then, we will compare the X-SFEM method with a $L^2$ projection method based on a deterministic X-FEM code (P-X-FEM). In example 2, we consider a plate with a random circular hole, randomness being modeled with three independent random variables. The X-SFEM method will be compared with Monte-Carlo method and P-X-FEM. In example 3, a simplified 2D modeling of a random welded joint will be considered. This last example presents a more complex random geometry whose representation by the level-set technique is obtained by suitable boolean operations on basic random domains. X-SFEM will be compared to a classical FEM approach and P-X-FEM.
## 7.1 Example 1: random plate in tension

### 7.1.1 Problem definition and approximation

<figure><img src="image_13.png" /><figcaption>Fig. 12. Example 1. Random plate in tension</figcaption></figure>

We consider the problem of a plate submitted to uniform tension, represented on figure 12. The plate lies in the random domain $\Omega(\theta) =]0, 1[\times] - \frac{a(\theta)}{2}, \frac{a(\theta)}{2}] \subset B =]0, 1[\times] - \frac{1}{2}, \frac{1}{2}]$, where $a(\theta)$ is a random variable writing: $a(\theta) = \xi_1(\theta) + \xi_2(\theta)$, where $\xi_1 \in U(0.75, 0.95)$ and $\xi_2 \in U(0, 0.05)$ are two statistically independent uniform random variables<fnref n="1" />. We will then work in the 2-dimensional stochastic domain $\Theta = [0.75, 0.95] \times [0, 0.05]$. We work under plane stress assumption and consider a homogeneous isotropic elastic material with Young modulus equal to 1 and Poisson coefficient equal to 0.3. The random domain is characterized by the following random level-set: for $x = (x, y) \in B$,

$$ \phi(\mathbf{x}, \theta) = \max\{y - \frac{a(\theta)}{2}, -\frac{a(\theta)}{2} - y\}. \qquad (81) $$

The plate is submitted to a uniform and unitary tension load $F = (\pm 1, 0)$ on $\Gamma_2(\theta) = \{\mathbf{x} \in \partial B; \phi(\mathbf{x}, \theta) < 0\}$. The exact solution $u$ to this problem writes: for $x = (x, y) \in \Omega(\theta)$, $u(x, \theta) = (x, -0.3y)$. Let us notice that $u$ does not depend on the elementary event $\theta$.

With the X-SFEM method, we work on a unique mesh $\mathcal{T}_h$ of $B$, represented on figure 13, where elements can be split into two groups: the first group ($e_i$) gathers all the elements surely in the domain and the second group ($e_c$) gathers all the elements possibly cut by the boundary. For the approximation at the stochastic level, we use a generalized polynomial chaos with degree $p = 1$, which should be sufficient regarding the exact solution.

In order to estimate the quality of the X-SFEM solution, we introduce the following global error indicator on the stress field:

$$ \varepsilon = \frac{\|\boldsymbol{\sigma}_{h,P} - \boldsymbol{\sigma}\|_{L^2(\mathcal{T}, d\mu)}}{\|\boldsymbol{\sigma}\|_{L^2(\mathcal{T}, d\mu)}}, \qquad (82) $$

<footnote n="1">¹ $\xi \in U(a, b)$ denotes a uniform random variable with values in $[a, b]$.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Example 1. X-SFEM mesh with 2 element groups: surely in the domain ($e_i$) and possibly cut by the boundary ($e_c$)</figcaption></figure>

where $\sigma$ and $\sigma_{h,P}$ respectively denote the exact and approximate stress field and where the $L^2$-norm is defined by

$$ \|\sigma\|_{L^2(\mathcal{P},d\mu)}^2 = \int_{\Theta} \int_{\Omega(\theta)} \sigma(\mathbf{x}, \theta) : \sigma(\mathbf{x}, \theta) \, dx \, dP(\theta). \quad (83) $$

We also introduce the error indicator $\varepsilon_K$, which is the local contribution of a finite element $K \in \mathcal{T}_h$ to the global error:

$$ \varepsilon_K = \frac{\|\sigma_{h,P} - \sigma\|_{L^2(\mathcal{P}_K, d\mu)}}{\|\sigma\|_{L^2(\mathcal{P}_K, d\mu)}}, \quad (84) $$

where $\mathcal{P}_K = \{(\mathbf{x}, \theta) \in B \times \Theta; \mathbf{x} \in K \cap \Omega(\theta)\}$ denotes the physical part of $K \times \Theta$.

### 7.1.2 Efficiency of the integration at stochastic level

First, we study the quality of the proposed integration technique at the stochastic level. In particular, we compare the two proposed recursive procedures for partitioning the stochastic domain $\Theta$ (see section 6.2.3). Let $K \in \mathcal{T}_h$ be a given finite element and let $C_K$ denote the associated cartesian mesh of $\Theta$, obtained by the isotropic or the anisotropic splitting procedure. We denote by $k$ the maximum order of the partition, which is directly related to the integration precision. Figure 14 presents stochastic domain partitions of order $k = 4$ which are associated with a particular finite element (shown on figure 12). This finite element belongs to the set of possibly cut elements ($e_c$). We observe that the number of stochastic subdomains is lower with the anisotropic splitting (69 versus 88). On this figure, the recursive splitting procedure has stopped for cells $C$ which are not colored in purple. We recall that the recursive split stops if all vertices $\xi \in \text{vertices}(C)$ correspond to the same state $S_K(\xi) = (\text{sign}(\phi_i(\xi)))_{i \in I_K}$ of the finite element $K$. If we increase the maximum order $k$ of the partition, only cells colored in purple will be split again.
Let us here recall that $\xi_1$ and $\xi_2$ have not the same range. Axis in Figure 14 have been rescaled for clarity.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Example 1. Stochastic partitions $\mathcal{C}_K$ of order $k = 4$ for a particular finite element $K$: isotropic splitting (left) and anisotropic splitting (right).</figcaption></figure>

Let $N_K$ be the number of subdomains in partition $\mathcal{C}_K$. Table 1 shows the mean value of $N_K$ over all $K \in (e_c)$, denoted by $\mu_N$, for both partition procedures and for different orders of partitions. We note that this number is lower with the anisotropic splitting. We denote by $T_K$ the computational time for creating partition $\mathcal{C}_K$ and integrating elementary quantities $\mathbf{A}_K$ and $\mathbf{b}_K$. Table 1 also shows the mean time $\mu_T$ over all elements $K \in (e_c)$. We note that this time is reduced with the anisotropic procedure, especially for high $k$ values. In fact, the computational time of X-SFEM mainly comes from the integration of element quantities for elements in $(e_c)$. The computational time $T_K$ for $K \in (e_i)$ and solving the Galerkin system of equations (56) is negligible for this problem (less than 1 second).

<table><caption>Table 1</caption><thead><tr><th>partition order</th><th></th><th>$k=0$</th><th>$k=1$</th><th>$k=2$</th><th>$k=3$</th><th>$k=4$</th><th>$k=5$</th></tr></thead><tbody><tr><td rowspan="2">$\mu_N$</td><td>isotropic</td><td>1</td><td>4</td><td>14</td><td>32</td><td>69</td><td>143</td></tr><tr><td>anisotropic</td><td>1</td><td>4</td><td>11</td><td>26</td><td>54</td><td>106</td></tr><tr><td rowspan="2">$\mu_T$</td><td>isotropic</td><td>0.04s</td><td>0.11s</td><td>0.24s</td><td>0.5s</td><td>1.1s</td><td>2.6s</td></tr><tr><td>anisotropic</td><td>0.04s</td><td>0.12s</td><td>0.22s</td><td>0.46s</td><td>1.0s</td><td>2.0s</td></tr></tbody></table>

Example 1. Comparison between the isotropic and anisotropic partitioning procedures for different orders $k$ of stochastic partitions: $\mu_N$ is the average number of stochastic subdomains for elements belonging to $(e_c)$ and $\mu_T$ is the corresponding average CPU times (for creating the stochastic partition and integrating element quantities).

Here, the approximation space contains the exact solution. Then, since X-SFEM is based on a Galerkin projection of the exact physical solution, it leads to the exact solution. However, the obtained approximate solution clearly depends on the partition order $k$ (related to the integration precision). Let us denote by $\boldsymbol{\sigma}_{h,P}^k$ the associate approximate stress field and by $\varepsilon^k$ the corresponding value of the global error indicator. On figure 15, we show the convergence
of $\epsilon^k$ with respect to $k$. We can observe a quick convergence towards the exact solution for both procedures. For instance, for $k = 2$, both procedures lead to a relative error inferior to $10^{-4}$.

<figure><img src="image_16.png" /><figcaption>Fig. 15. Example 1. Integration precision: convergence of global error indicator $\epsilon^k$ with respect to $k$ for isotropic and anisotropic procedures.</figcaption></figure>

### 7.1.3 Quality of the X-SFEM solution

Some examples of post-processing will now illustrate the quality of the solution obtained by X-SFEM. The solution being explicit in terms of basic random variables, post-processing can be performed at a very low cost. The following results have been obtained with the mesh presented on figure 13 and with a stochastic integration based on a stochastic partition of order $k = 2$. Figure 16 presents the response surfaces of the horizontal displacement for two points: $\boldsymbol{x}_1 = (1, 0.5)$ which is surely inside the random domain and $\boldsymbol{x}_2 = (1, 0.95)$ which is possibly inside or outside the random domain. For both points, we observe that the approximation matches very well the exact physical solution. For point $\boldsymbol{x}_2$, we observe that the non-physical part of the solution is the unique possible prolongation (in the approximation space) of the physical part of the solution.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Example 1. Response surfaces of horizontal displacement for points $\boldsymbol{x}_1 = (1, 0.5)$ and $\boldsymbol{x}_2 = (1, 0.95)$.</figcaption></figure>

Figure 17 represents the stress field component $\sigma_{xx}$ for a particular event such as $a(\theta) = 0.75$. The iso-zero of the corresponding outcome $\phi(\boldsymbol{x}, \theta)$ of the level-set is represented by the red line. We observe that the stress solution is
uniform and equal to 1 according to the exact solution even in the non-physical
domain. We also show on the same figure the local contribution $\epsilon_K$ of each
finite element to the global error. We note that these local contributions are
really low and concentrated on few elements.

<figure><img src="image_18.png" /><figcaption>Fig. 17. Example 1. Stresses $\sigma_{xx}$ (left) and local contributions to the error $\epsilon_K^k$ for $k=2$ (right) for a particular event ($a(\theta) = 0.75$).</figcaption></figure>

### 7.1.4 Comparison with a projection method based on X-FEM (P-X-FEM)

We now focus on the comparison between the X-SFEM method and P-X-FEM
(see section 5.6.1). Here, let us denote by $\bar{\sigma}_{h,P}$ the solution obtained by P-X-
FEM and by $\tilde{\epsilon}^p$ and $\tilde{\epsilon}_K^p$ the corresponding values of the global and local error
indicators (defined in (82) and (84)). Superscript $p$ is related to the chosen
degree of the generalized polynomial chaos. Figure 18 shows the evolution of
$\tilde{\epsilon}^p$ with respect to $p$. We observe a very slow convergence of P-X-FEM. Even
for a degree $p=8$, the error is still high as opposed to the X-SFEM solution
which reaches a very high precision with a degree equal to 1. This is due to
the bad definition of the P-X-FEM approximation, which gives weight to the
non-physical part of the solution.

<figure><img src="image_19.png" /><figcaption>Fig. 18. Example 1. Convergence of P-X-FEM: error indicator $\tilde{\epsilon}^p$ with respect to generalized polynomial chaos degree $p$.</figcaption></figure>

As we can see on figure 19, which plots the local error indicator, the error
is concentrated on elements which are possibly cut by the boundary, more
precisely on elements which have a non-zero probability to be outside the
random domain. Figure 19 also plots the stress field component $\sigma_{xx}$ for a
particular event such that $a(\theta) = 0.75$.
<figure><img src="image_20.png" /><figcaption>Fig. 19. Example 1. P-X-FEM: stress field component $\sigma_{xx}$ (left) and local contribution $\tilde{\varepsilon}_K^p$ to the global error (right) for a particular event ($a(\theta) = 0.75$)</figcaption></figure>

Figure 20 shows the horizontal displacement for points $\boldsymbol{x}_1$ and $\boldsymbol{x}_2$ for which the X-SFEM solutions were presented on figure 16. Those results have been calculated with a generalized polynomial chaos of degree $p = 8$. We observe that the P-X-FEM solution for $\boldsymbol{x}_1$ (point surely inside the domain) is very good, as the X-SFEM solution. However it is not the case for point $\boldsymbol{x}_2$, which is possibly outside the domain: that clearly illustrates the bad behavior of the approximation based on a $L^2$ projection.

<figure><img src="image_21.png" /><figcaption>Fig. 20. Example 1. P-X-FEM: response surfaces of horizontal displacements for points $\boldsymbol{x}_1$ and $\boldsymbol{x}_2$.</figcaption></figure>

### 7.1.5 A particular case where no specific integration is required

At the spatial level, if we use a coarse finite element mesh as presented on figure 21, we can observe that the X-SFEM leads the exact solution in the numerical sense (relative error less than $10^{-14}$). This can be explained as follows. For each finite element $K$, the $S_K(\xi) = (sign(\phi_i(\xi)))_{i \in I_K}$ is the same for all elementary events $\xi \in \Theta$ (illustrated on figure 21, which plots all possible outcomes of the random level-set). Therefore, local quantities $A_K(\xi)$ and $b_K(\xi)$ are global polynomial (of degree 2 in $\xi$) on the whole stochastic domain (see section 6.2.2). In this case, a classical quadrature rule defined on the whole stochastic domain $\Theta$ (without partitioning) leads to an exact integration of element quantities. Then, we obtain the exact approximate solution of X-SFEM, which is the exact solution for this simple problem.
<figure><img src="image_22.png" /><figcaption>Fig. 21. Example 1. Coarse X-SFEM finite element mesh.</figcaption></figure>

## 7.2 Example 2: plate with a random circular hole in tension

We consider the problem of a plate with a random circular hole submitted to uniaxial tension. We consider a homogeneous isotropic elastic material, under plane stress assumption, with Young modulus equal to 1 and Poisson coefficient equal to zero. The plate lies in a square domain $B = ]0, 3[ \times ]0, 3[$. The location of the center of the hole $c = (c_x, c_y)$ depends on two statistically independent uniform random variables $c_x(\theta) \in U(1.4, 1.5)$ and $c_y(\theta) \in U(1.5, 1.55)$. The radius is also modeled with a uniform random variable $R(\theta) \in U(0.45, 0.5)$, which is statistically independent of $c_x$ and $c_y$. In this case, the level-set function can be explicitly written as follows: $\phi(\mathbf{x}, \theta) = R(\theta) - \|\mathbf{x} - \mathbf{c}(\theta)\|$. With the X-SFEM method, we use a unique mesh represented on figure 22. The plate is submitted to a uniform deterministic tension load: $\mathbf{F} = (-1, 0)$ on $\{\mathbf{x} \in \partial B; x = 0\}$ and $\mathbf{F} = (1, 0)$ on $\{\mathbf{x} \in \partial B; x = 3\}$.

<figure><img src="image_23.png" /><figcaption>Fig. 22. Example 2. Plate with a random hole in tension; X-SFEM mesh with three element groups: surely inside the domain ($e_i$), possibly cut by the boundary of the hole ($e_c$) and surely outside the domain ($e_o$).</figcaption></figure>

We denote by $\boldsymbol{\xi} = (c_x, c_y, R)$. The stochastic domain is then $\boldsymbol{\Theta} = [1.4, 1.5] \times [1.5, 1.55] \times [0.45, 0.5]$. At the stochastic level, we use for the approximation a 3-dimensional generalized polynomial chaos with degree $p = 2$.
### 7.2.1 Efficiency of the integration at stochastic level

We study the efficiency of the stochastic integration technique based on a partitioning of the stochastic domain (see section 6.2.3). Figure 23 presents different partitions $C_K$ obtained with the isotropic and anisotropic procedures for a particular finite element $K$ (shown on figure 22) and for different orders $k$ of the partitions. We recall that cells which are not colored in purple correspond to cells whose vertices $\xi \in \Theta$ have all the same state $S_K(\xi)$. Then, for these cells, the splitting procedure has reached its stopping criterium.

Let $N_K$ be the number of subdomains in partition $C_K$. Table 2 indicates the mean value of $N_K$ over all $K \in (e_c)$, denoted by $\mu_N$, for both partition procedures and for different orders of partitions. It also indicated the mean time $\mu_T$, over all elements $K \in (e_c)$, which is required for creating stochastic partitions and integrating element quantities $A_K$ and $b_K$. As in the previous example, we notice that $\mu_N$ is reduced with the anisotropic procedure: for instance, with $k=4$ or $k=5$, this number is up to 20% lower than with the isotropic splitting. Therefore, it also leads to a reduction of computation time $\mu_T$.

<figure><img src="image_24.png" /><figcaption>Fig. 23. Example 2. Stochastic partitions $C_K$ associated with a finite element $K$ for different orders $k$: comparison between isotropic splitting (left) and anisotropic splitting (right).</figcaption></figure>

Let us denote by $\varepsilon^k$ the global $L^2$ error indicator between the approximate stress field $\sigma_{h,P}^k$ (for a given $k$) and a reference approximate stress field $\sigma_{h,P}$
<table><caption>Table 2</caption><thead><tr><th colspan="2">partition order</th><th>k = 1</th><th>k = 2</th><th>k = 3</th><th>k = 4</th><th>k = 5</th></tr></thead><tbody><tr><td rowspan="2">μ<sub>N</sub></td><td>isotropic</td><td>8</td><td>50</td><td>245</td><td>1076</td><td>4493</td></tr><tr><td>anisotropic</td><td>7</td><td>41</td><td>198</td><td>840</td><td>3413</td></tr><tr><td rowspan="2">μ<sub>T</sub></td><td>isotropic</td><td>0.45s</td><td>2.3s</td><td>12.1s</td><td>49.8s</td><td>210s</td></tr><tr><td>anisotropic</td><td>0.43s</td><td>2.1s</td><td>10.3s</td><td>42.4s</td><td>188s</td></tr></tbody></table>

Example 2. Comparison between the isotropic and anisotropic partitioning procedures for different orders $k$ of stochastic partitions: $\mu_N$ is the average number of stochastic subdomains for elements belonging to $(e_c)$ and $\mu_T$ is the corresponding average CPU times (for creating the stochastic partition and integrating element quantities).

(computed with an integration on an isotropic partition of order $k = 6$). On figure 24, we can observe a good convergence of $\varepsilon^k$ with respect to $k$ for both partition procedures. In the following, we will use the X-SFEM approximate solution corresponding to an anisotropic splitting of order $k = 3$, which leads to a relative error lower than $10^{-3}$ on the calculation of the approximate solution.

<figure><img src="image_25.png" /><figcaption>Fig. 24. Example 2. Convergence of $\varepsilon^k$ with respect to $k$ for the isotropic and anisotropic partitions of $\Theta$.</figcaption></figure>

### 7.2.2 Comparison with classical approaches

In this section, we compare the X-SFEM solution with solutions calculated with classical approaches: figure 25 shows a comparison between FEM and X-SFEM solutions for a particular outcome of the geometry. We observe a very good matching between stress fields $\sigma_{xx}$.

Figure 26 shows the probability density function for the horizontal elongation $\Delta U(\theta) = u_x(\mathbf{x}_2, \theta) - u_x(\mathbf{x}_1, \theta)$ where $\mathbf{x}_1 = (0, 1.5)$ and $\mathbf{x}_2 = (3, 1.5)$. The reference solution has been calculated with a Monte-Carlo approach with $10^4$ samplings using an X-FEM code to solve each deterministic problem. We observe a good X-SFEM solution according to the reference solution, even
<figure><img src="image_26.png" /><figcaption>Fig. 25. Example 2. Comparison between FEM (left) and X-SFEM (right): stresses $\sigma_{xx}$ for $c_x = 1.4$, $c_y = 1.5$ and $R = 0.45$. (for vizualisation, stresses in elements outside the physical domain are set to zero)</figcaption></figure>

with a polynomial chaos of degree $p = 2$. Finally, let us indicate that on a personal computer, the computational time took about 180 minutes with the Monte-Carlo approach and only about 6 minutes with the X-SFEM method.

<figure><img src="image_27.png" /><figcaption>Fig. 26. Example 2. Probability density function of $\Delta U$: comparison of X-SFEM solutions (for different polynomial chaos degree $p$) with a Monte-Carlo approach coupled with X-FEM ($10^4$ samplings).</figcaption></figure>

### 7.3 Example 3: random welded joint

We consider the problem of a structure with a random welded joint. The problem is represented on figure 27. The structure is embedded on a part $\Gamma_1$ of its boundary and submitted a uniform bending surface load $\mathbf{F} = (0, 0.1) \, \text{MPa}$ on $\Gamma_2$. The shape of the welded joint is represented by an elliptical arch characterized by its two semi-axes $a$ and $b$ which also define the center of the ellipse. Those two parameters are taken as statistically independent uniform random variables $a(\theta) \in U(0.05, 0.15)$ and $b(\theta) \in U(0.05, 0.15)$. We then work in a 2-dimensional stochastic domain $\mathbf{\Theta} = [0.05, 0.15] \times [0.05, 0.15]$. We consider that the material is homogeneous (same material for the structure and
<figure><img src="image_28.png" /><figcaption>Fig. 27. Example 3. Welded joint problem</figcaption></figure>

the welded part) isotropic elastic with Young modulus equal to 210 GPa and Poisson coefficient equal to 0.3. We solve the problem under the plane strain assumption. With the X-SFEM method, we use a single mesh lying in the domain $B = ]0, 0.4[ \times ]0, 0.7[$. In order to represent the random geometry, we introduce boolean operations on several simple random domain (ellipses and semi-planes), characterized by basic level-sets (see section 5.2.2). Figure 28 presents the X-SFEM finite element mesh composed by 2304 3-nodes triangles with average size $h = 1.5 \ 10^{-2}$. Figure 29 shows the interpolation of the level-set for an outcome of the geometry corresponding to $a = b = 0.1$.

<figure><img src="image_29.png" /><figcaption>Fig. 28. Example 3. X-SFEM mesh with three element groups: surely inside the domain ($e_i$), possibly cut by the boundary ($e_c$) and surely outside the domain ($e_o$).</figcaption></figure>

### 7.3.1 Efficiency of the integration at stochastic level

As in previous examples, we first study the influence of the integration technique at the stochastic level. At the stochastic level, we use a 2-dimensional
<figure><img src="image_30.png" /><figcaption>Fig. 29. Example 3. Outcome of the random level-set corresponding to $a = 0.1$ and $b = 0.1$.</figcaption></figure>

generalized polynomial chaos of degree $p = 2$. Let $k$ denote the maximum order of stochastic partitions (anisotropic or isotropic) and let us introduce the global $L^2$ error indicator $\varepsilon^k$ between the approximate stress field $\sigma_{h,P}^k$ (for a given order $k$ of stochastic partitions) and a reference approximate stress field $\sigma_{h,P}$ (computed with an integration on an anisotropic partition of order $k = 6$). This error indicator is defined as in (82). Let us also introduce the local error contribution $\varepsilon_K^k$ of a finite element $K$, defined as in (84). Figure 30 shows the convergence of $\varepsilon^k$ with respect to $k$ and figure 31 shows the local errors $\varepsilon_K^k$ for $k = 3$ and $k = 4$. We note a good convergence towards the reference solution and, as we have noticed in previous examples, an order $k$ equal to 3 seems to be sufficient to perform a good integration at the stochastic level. We still observe that errors are concentrated on a few elements which are possibly cut by the boundary.

<figure><img src="image_31.png" /><figcaption>Fig. 30. Example 3. Convergence of $\varepsilon^k$ with respect to the order $k$ of stochastic partitions (isotropic and anisotropic).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 31. Example 3. Local errors $\epsilon_K^k$ for stochastic partitions orders $k = 3$ and $k = 4$.</figcaption></figure>

### 7.3.2 Comparison between X-SFEM, P-X-FEM and a classical FEM approach

In this section, we focus on comparison between the X-SFEM solution and the solutions obtained with two other approaches: a classical FEM approach with remeshings and P-X-FEM ($L^2$ projection method based on X-FEM). X-SFEM and P-X-FEM use the same finite element mesh but we use a generalized polynomial chaos of degree $p = 4$ for X-SFEM and of degree $p = 8$ for P-X-FEM. FEM solutions have been obtained with conforming meshes, corresponding to given outcomes of the geometry, with a mesh size equivalent to the X-SFEM and P-X-FEM mesh size. First we compare the solutions for the vertical displacement at two points $P_1$ and $P_2$, shown on figure 27, which are respectively surely inside the domain and possibly inside the domain. Figures 32 and 33 show the response surfaces corresponding to each approach for those two points. For point $P_1$, we observe a very good agreement between the response surfaces of the three approaches. For point $P_2$, X-SFEM and P-X-FEM solutions are composed by a physical part and a non-physical part. We observe that X-SFEM (based on a Galerkin projection of the physical part of the exact solution) is much better than P-X-FEM (based on a $L^2$ projection of an artificial prolongation of the exact solution).

<figure><img src="image_33.png" /><figcaption>Fig. 32. Example 3. Response surfaces for vertical displacement at point $P_1$ surely inside the domain: comparison between FEM (left), X-SFEM (center) and P-X-FEM (right).</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Fig. 33. Example 3. Response surfaces for vertical displacement at point $P_2$ possibly inside the domain: comparison between FEM (left), X-SFEM (center) and P-X-FEM (right).</figcaption></figure>

Figure 34 shows the response surfaces of the maximum value of stress field component $\sigma_{xx}$ for the three approaches. This quantity of interest, writing $\sigma_{xx}^{max}(\xi) = \max_{x \in \Omega(\xi)}(\sigma_{xx}(x, \xi))$, can be very useful in practice. On this figure, we still observe a good agreement between FEM and X-SFEM response surfaces. P-X-FEM solution is clearly less accurate than X-SFEM. While FEM and X-SFEM give a maximum stress for $a = 0.05$ and $b = 0.15$, P-X-FEM gives this maximum at $a = 0.15$ and $b = 0.15$. Moreover, the maximum stress value $sup_{\xi \in \Theta}(\sigma_{xx}^{max}(\xi))$ given by P-X-FEM is much higher than the one given by FEM and X-SFEM. Finally, let us indicate the computational times required by these approaches on a personal computer: 5 minutes for X-SFEM, 13 minutes for P-X-FEM and 4 hours for a Monte-Carlo approach coupled with a classical FEM code and $10^4$ samplings.

<figure><img src="image_35.png" /><figcaption>Fig. 34. Example 3. Response surfaces for maximum stress $\sigma_{xx}^{max}$: comparison between FEM (left), X-SFEM (center) and P-X-FEM (right).</figcaption></figure>

### 7.3.3 Convergence of the X-SFEM solution

In this section, we focus on the convergence of the X-SFEM solution. We introduce the following global error indicator between the X-SFEM displacement solution $u_{h,P}$ and a reference displacement field $u$:

$$ \varepsilon_{h,p} = \frac{\|u_{h,P} - u\|_{L^2(\mathcal{P},d\mu)}}{\|u\|_{L^2(\mathcal{P},d\mu)}}, \qquad (85) $$

where index $p$ denotes the degree of the generalized polynomial chaos, $h$ the average finite element mesh size and $\| \cdot \|_{L^2(\mathcal{P},d\mu)}$ the $L^2$ norm on the physical
part of the displacement field, defined as for stresses in equation (83). The reference solution **u** is obtained with X-SFEM with a mesh size *h* ≈ 0.015 and a degree *p* = 6. Figure 35 shows the convergence of ε<sub>h,p</sub> with respect to *p* and *h*. For a given *h*, we observe a very quick convergence with respect to *p*. For *p* ≥ 2, we can conclude that the contribution to error comes essentially from the finite element error. For a given *p*, we observe a classical convergence in O(*h*<sup>2</sup>).

<figure><img src="image_36.png" /><figcaption>Fig. 35. Example 3. Convergence of the X-SFEM approximation: error indicator εh,p with respect to p (top) and h (bottom)</figcaption></figure>

8 Conclusion

We proposed a Stochastic Finite Element Method to solve stochastic par-
tial differential equations defined on random domains. The first point of this
method lies in the use of the level-set technique, which represents in an im-
plicit way the geometry and allow the handling of complex random geometries.
The second point lies in the construction of a tensor product approximation
space, which is made possible by considering prolongation of solutions on a
deterministic spatial domain. A Galerkin criterium is then used to define the
approximation. Another definition of the approximation has been also intro-
duced, which is based on a classical L² projection method based on a classical
X-FEM code. The superiority of the Galerkin projection has been illustrated
and interpreted from a mathematical point of view. The technical aspects for
the construction and the resolution of the discretized problem have been de-
tailed. A particular care has been devoted to the numerical integration of the weak form, requiring the development of a specific quadrature technique at the stochastic level. For the numerical examples treated in this article, the proposed stochastic quadrature technique has given very good results since elementary matrices and vectors had a nice dependence on the random variables. The case where these quantities have a more complex dependence on the random variables and also the case of higher stochastic dimension are key questions which are currently under investigation. These points will certainly require the development of simpler quadrature techniques and the use of the high degree of parallelism of the method in order to reduce computational costs.

## Acknowledgement

This work is supported by the French National Research Agency (grant ANR-06-JCJC-0064) and by the European Community and FEDER funds within the project MEDACHS Interreg IIIB.

## References

[1] R. G. Ghanem and P. D. Spanos. *Stochastic finite elements: a spectral approach*. Springer, Berlin, 1991.

[2] I. Babuška, R. Tempone, and G. E. Zouraris. Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation. *Computer Methods in Applied Mechanics and Engineering*, 194:1251–1294, 2005.

[3] H. G. Matthies and A. Keese. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 194(12-16):1295–1331, 2005.

[4] R. G. Ghanem and R. M. Kruger. Numerical solution of spectral stochastic finite element systems. *Computer Methods in Applied Mechanics and Engineering*, 129:289–303, 1996.

[5] R. G. Ghanem and W. Brzakala, Stochastic finite-element analysis of soil layers with random interface. *Journal of Engineering Mechanics*, 122(4):361–369, 1996

[6] M. F. Pellissetti and R. G. Ghanem. Iterative solution of systems of linear equations arising in the context of stochastic finite elements. *Advances in Engineering Software*, 31:607–616, 2000.
[7] A. Keese and H. G. Mathhies. Hierarchical parallelisation for the solution of stochastic finite element equations. *Computer Methods in Applied Mechanics and Engineering*, 83:1033-1047, 2005.

[8] A. Nouy. A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations. *Computer Methods in Applied Methods in Engineering*, 196(45-48):4521-4537, 2007.

[9] R. E. Caflisch. Monte carlo and quasi-monte carlo methods. *Acta. Numer.*, 7:1-49, 1998.

[10] M. Papadrakakis and V. Papadopoulos. Robust and efficient methods for stochastic finite element analysis using monte carlo simulation. *Computer Methods in Applied Mechanics and Engineering*, 134:325-340, 1996.

[11] B. Puig, F. Poirion, and C. Soize. Non-gaussian simulation using hermite polynomial expansion: convergences. *Probabilistic Engineering Mechanics*, 17:253-264, 2002.

[12] M. Berveiller, B. Sudret, and M. Lemaire. Stochastic finite element: a non intrusive approach by regression. *European Journal of Computational Mechanics*, 15:81-92, 2006.

[13] D. B. Xiu and D. M. Tartakovskiy. Numerical methods for differential equations in random domains. *SIAM J. Sci. Comput.*, 28(3):1167-1185, 2006.

[14] D. M. Tartakovskiy and D. B. Xiu. Stochastic analysis of transport in tubes with rough walls. *Journal of Computational Physics*, 217:248-259, 2006.

[15] A. Nouy, F. Schoefs, and N. Moës. X-SFEM, a computational technique based on X-FEM to deal with random shapes. *European Journal of Computational Mechanics*, 16(2):277-293, 2007.

[16] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. *Int. J. for Numerical Methods in Engineering*, 46:131-150, 1999.

[17] T. Belytschko, C. Parimi, N. Moës, N. Sukumar, and S. Usui. Structured extended finite element methods for solids defined by implicit surfaces. *Int. J. for Numerical Methods in Engineering*, 56:609-635, 2003.

[18] J.A. Sethian. *Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science*. Cambridge University Press, Cambridge, UK, 1999.

[19] C. Daux, N. Moës, J. Dolbow, N. Sukumar, and T. Belytschko. Arbitrary branched and intersecting cracks with the extended finite element method. *Int. J. for Numerical Methods in Engineering*, 48:1741-1760, 2000.

[20] N. Sukumar, D.L. Chopp, N. Moës, and T. Belytschko. Modeling holes and inclusions by level sets in the extended finite-element method. *Computer Methods in Applied Mechanics and Engineering*, 190:6183-6200, 2001.
[21] J. M. Melenk and I. Babuška. The partition of unity method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering*, 39:289–314, 1996.

[22] I. Ramière, P. Angot, and M. Belliard. A fictitious domain approach with spread interface for elliptic problems with general boundary conditions. *Computer Methods in Applied Mechanics and Engineering*, 196(4-6):766–781, 2007.

[23] R. Glowinski, T. Pan, and J. Périaux. A fictitious domain method for dirichlet problems and applications. *Computer Methods in Applied Mechanics and Engineering*, 111:283–303, 1994.

[24] V. Girault and R. Glowinski. Error analysis of a fictitious domain method applied to a dirichlet problem. *Japan J. Indust. Appl. Math.*, 12:487–514, 1995.

[25] R. Glowinski and Y. Kuznetsov. On the solution of the dirichlet problem for linear elliptic operators by a distributed lagrange multiplier method. *C. R. Acad. Sci. Paris, Série I*, 327:693–698, 1998.

[26] C. Canuto and T. Kozubek. A fictitious domain approach to the numerical solution of pdes in stochastic domains. *Numerische Mathematik*, 107(2):257–293, 2007.

[27] C. Soize. Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators. *Computer Methods in Applied Mechanics and Engineering*, 195(1-3):26–64, 2006.

[28] I. Babuška and J. Chleboun. Effects of uncertainties in the domain on the solution of neumann boundary value problems in two spatial dimensions. *Mathematics of Computation*, 71(240):1339–1370, 2002.

[29] I. Babuška and J. Chleboun. Effects of uncertainties in the domain on the solution of dirichlet boundary value problems. *Numer. Math.*, 93(4):583–610, 2003.

[30] M. Loève. *Probability Theory. I, fourth edition, in: Graduate Texts in Mathematics*, vol. 45. Springer-Verlag, New York, 1977.

[31] M. Loève. *Probability Theory. II, fourth edition, in: Graduate Texts in Mathematics*, vol. 46. Springer-Verlag, New York, 1978.

[32] B. Øksendal. *Stochastic Differential Equations. An Introduction with Applications, fifth ed.* Springer-Verlag, 1998.

[33] P. Frauenfelder, C. Schwab, and R. A. Todor. Finite elements for elliptic problems with stochastic coefficients. *Computer Methods in Applied Mechanics and Engineering*, 194(2-5):205–228, 2005.

[34] N. Wiener. The homogeneous chaos. *Am. J. Math.*, 60:897–936, 1938.

[35] D. B. Xiu and G. E. Karniadakis. The Wiener-Askey polynomial chaos for stochastic differential equations. *SIAM J. Sci. Comput.*, 24(2):619–644, 2002.
[36] C. Soize and R. Ghanem. Physical systems with random uncertainties: chaos representations with arbitrary probability measure. *SIAM J. Sci. Comput.*, 26(2):395–410, 2004.

[37] M. Deb, I. Babuška, and J. T. Oden. Solution of stochastic partial differential equations using galerkin finite element techniques. *Computer Methods in Applied Mechanics and Engineering*, 190:6359–6372, 2001.

[38] O. P. Le Maître, O. M. Knio, H. N. Najm, and R. G. Ghanem. Uncertainty propagation using Wiener-Haar expansions. *Journal of Computational Physics*, 197(1):28–57, 2004.

[39] O. P. Le Maître, H. N. Najm, R. G. Ghanem, and O. M. Knio. Multi-resolution analysis of wiener-type uncertainty propagation schemes. *Journal of Computational Physics*, 197(2):502–531, 2004.

[40] J. Thompson, Z. Warsi, and C. Mastin. *Numerical Grid Generation*. Elsevier Science, Amsterdam, 1985.

[41] F. Brezzi. On the existence, uniqueness and approximation of saddlepoint problems arising from lagrange multipliers. *RAIRO Anal. Numér.*, R-2:129–151, 1974.

[42] N. Moës, E. Béchet, and M. Tourbier. Imposing dirichlet boundary conditions in the extended finite element method. *Int. J. for Numerical Methods in Engineering*, 67(12):1641–1669, 2006.

[43] S. Géniaut, P. Massin, and N. Moës. A stable 3d contact formulation using x-fem. *European Journal of Computational Mechanics*, 16(2):259–275, 2007.

[44] I. Babuška. The finite element method with penalty. *Mathematics of Computation*, 27(122):221–228, 1973.

[45] J. A. Nitsche. Uber ein variationsprinzip zur Lösung dirichlet-problemen bei verwendung von Teilräumen, die keinen randbedingungen unteworfen sind. *Abh. Math. Sem. Univ. Hamburg*, 36:9–15, 1971.

[46] R. G. Ghanem. Ingredients for a general purpose stochastic finite elements implementation. *Computer Methods in Applied Mechanics and Engineering*, 168:19–34, 1999.