Numerical upscaling of parametric microstructures in a possibilistic uncertainty framework with tensor trains

A fuzzy arithmetic framework for the efficient possibilistic propagation of shape uncertainties based on a novel fuzzy edge detection method is introduced. The shape uncertainties stem from a blurred image that encodes the distribution of two phases in a composite material. The proposed framework employs computational homogenisation to upscale the shape uncertainty to a effective material with fuzzy material properties. For this, many samples of a linear elasticity problem have to be computed, which is significantly sped up by a highly accurate low-rank tensor surrogate. To ensure the continuity of the underlying mapping from shape parametrisation to the upscaled material behaviour, a diffeomorphism is constructed by generating an appropriate family of meshes via transformation of a reference mesh. The shape uncertainty is then propagated to measure the distance of the upscaled material to the isotropic and orthotropic material class. Finally, the fuzzy effective material is used to compute bounds for the average displacement of a non-homogenized material with uncertain star-shaped inclusion shapes.


Introduction
Composite materials are ubiquitous in many applications. Whether they are formed by chance (like with unfavourable impurities in metal) or by design (like pebbles in concrete), the mechanical properties of the whole material are determined by the resulting composite structure [1]. The different Robert  Institut für Geometrie und Praktische Mathematik, RWTH Aachen, Templergraben 55, Aachen 52056, Germany types of composite materials-e.g. metal matrix composites [2], fiber-reinforced polymers [3], composite wood [4] and other advanced composite materials [5]-find their use in a wide range of applications, e.g. masonry, aerospace industry, wind power plants and sports equipment. In many applications, the added composite material improves a base matrix material in terms of wear resistance, damping properties and mechanical strength, while keeping the same weight. However, such improved composite materials are a result of empirical studies, experimentation and chance. It hence is obvious that a thorough understanding of composite materials via theoretical and accurate computational models is desirable to predict and systematically improve the properties and applicability of composites.
What makes this task challenging are the influence of material behaviour on the micro-and macro scales, as well as the multifaceted sources of uncertainties. For instance, the involved length scales may span up to ten orders of magnitudes, i.e. the size of a nano particle of 10 −6 m is embedded in a material with a length scale of 10 −1 −10 1 m [2]. The same holds for the material constants if the constituents are fundamentally different as in metals and polymers. A standard finite element approach becomes very costly in such a setting since the details of the composite have to be resolved by the mesh for accurate simulations. Additionally, one has to handle the uncertainty of the material constants such as the position, form and size of the inclusions. This increases the costs even further since many realisations are required to correctly capture the uncertainty statistically. Each new realisations requires a costly (possible automatic) remeshing of the computational domain. The mapping from geometric parameters-namely position and shape of each inclusiononto the generated meshes are usually discontinuous. This introduces unwanted effects that we tackle in Sect. 4.1.
The first challenge-the span of scales-is tackled by homogenisation methods, a term introduced by Babuška in [6]. In some form or another, these methods incorporate micro-scale behaviour into an adjusted macro-scale model, replacing the full composite material model by a corrected homogeneous one, obtained by an asymptotic limit of an assumed (periodic) domain. The central idea is to derive equations that describe the effective material properties analytically. Alternatively, computational approaches have been devised to solve particular micro-scale problems to deduce the adjusted macro-scale behaviour numerically [7,8]. Some notable examples are stochastic homogenisation, see [9][10][11][12], projection based homogenisation [13] and (stochastic) representative volume element methods [14].
The second challenge-the uncertainty of the material-is sometimes neglected [15,16] by only considering deterministic material properties. If the probability distribution of the material constants is known, the uncertainty can be modelled with precise probabilities [17]. Ignoring the uncertainty seems valid if it has little influence on the system's response or the model data is known sufficiently accurate and is free of inherent fluctuations. Using precise probabilities is valid if the distributions of the material constants are known precisely. However, as was pointed out by Motamed and Babuška [18], stochastic models based on precise probabilities are not always able to model the uncertainty in composite materials. Instead, they propose a model based on an imprecise probability theory. Examples for this are evidence theory [19], random set theory [20], possibility theory [21] andmore recently-optimal Uncertainty Quantification (UQ) [22]. In comparison to precise probability, imprecise probability methods are able to provide estimates of uncertainties based only on a small set of data and few assumptions, which would be far too limited for a probabilistic method. For an illuminating work that dissects the difference between precise and imprecise probabilities, we refer the reader to [23]. As a motivation the approach developed in what follows, we especially point to the false confidence theorem, which qualitatively states that "probability dilution is a symptom of a fundamental deficiency in probabilistic representations of statistical inference, in which there are propositions that will consistently be assigned a high degree of belief, regardless of whether or not they are true" [23].
In this work, we consider a practical problem that is supposed to illustrate when imprecise probabilities are appropriate (in particular more so than probabilistic methods), the challenges performing a non-probabilistic uncertainty propagation and quantification, and ways to overcome the computational difficulties. We assume a material consisting of two phases, namely a soft matrix phase and a hard inclusion phase. Both phases are linear elastic with precisely known Young's modulus and Poisson ratio. The inclusions repeat periodically in a checkerboard fashion. Consequently, a classical numerical homogenisation method would yield a homogenised (globally constant) macro-scale material. Considering real applications, the main problem stems from the unknown shape of the inclusions. Often, the shape is retrieved with a non-intrusive imaging process, e.g. computed tomography scans or magnetic resonance imaging. To illustrate this, Fig. 1 shows a computed tomography scan of an imperfect adhesive bond within a Henkel beam, cf. [24]. Apparently, the image is noisy, blurred, pixelated and exhibits artefacts. While it is possible to de-noise and de-blurr the image with image post-processing, which might make the identification of inclusions possible, different numerical methods yield different shapes [25]. This reveals the shape uncertainty inherent to the image. Here, we just focus on one blurred image per inclusion, assuming that the artefacts and noise were successfully removed already, see Fig. 3 for an example. Opposite to this, in the case of multiple CT scans of the same inclusion, the additional information can be employed to reduce the shape uncertainty and thus refine the epistemic uncertainty model until the uncertainty becomes irreducible and hereby aleatoric uncertainty models are preferable.
To tackle a problem in this setting, we introduce a novel fuzzy edge detection based on possibility theory, presented in Sect. 2.1 and a restriction of possible interface boundaries via bounded total variation in Sect. 3. This yields a fuzzy model of the boundary which in turn introduces a computational challenge in terms of a numerical homogenisation problem as discussed in Sect. 2.2. For this, the micro-model has to be solved very often in order to propagate the uncertainty of the boundary to the homogenised material model. To alleviate this expensive task, we introduce a highly accurate rank adaptive low-rank tensor surrogate in Sect. 4. In the numerical experiments Sect. 5, the surrogate model is validated numerically. Moreover, we measure the distance of the homogenised material to the class of isotropic and orthotropic elastic tensors and eventually use the homogenised material to perform a worst/best case analysis for a full matrix composite model with 64 inclusions, i.e., with help of the homogenised material we will find bounds for the average displacement of the non-homogenised material. Fig. 1 Computed tomography scan of a Henkel beam taken from [24]. Material imperfections are noticeable but the present resolution does not reveal the exact shapes of the imperfections

Basics
This section serves as a brief introduction of three fundamental topics that form the basis of this work. First, fuzzy set theory which is used to model the uncertainty ingrained in the blurred image. Second, the computational up-scaling method (numerical homogenisation). Third, classes of constitutive tensors to measure the distance of the up-scaled material to the isotropic and orthotropic material are discussed.
Fuzzy sets are common sets that are equipped with a membership function which assigns each element in the set a value between zero and one. This values only task is to communicate a degree of belongingness to the set. The meaning of this value depends on the problem and-more importantlyon the community using this value to formalise uncertainty. Arguably, this ambiguity is a feature and not a short-coming of this theory since it avoids assumptions that are usually made by other approaches. In probability theory for instance, one has to assume a prior distribution before updating the posterior distribution with new samples. The following definition introduces the terminology of fuzzy sets.

Definition 1 (Fuzzy set/variable, α-cuts and interactivity)
Let Z = ∅ be a set and μ : Z → [0, 1] be a map such that there exists z ∈ Z with μ(z) = 1. The map μ is called (normalised) membership function. We define a (normalised) fuzzy setz on Z bỹ If μ(Z ) = {0, 1} with unique z * ∈ Z and μ(z * ) = 1 thenz is called a crisp set. Furthermore, we denote by F(Z ) the set of all fuzzy sets on Z . We thus simply writez ∈ F(Z ). For each fuzzy setz ∈ F(Z ) we denote the associated membership function by μz. If Z ⊂ R N for N ∈ N, we callz ∈ F(Z ) a fuzzy variable (N = 1) or vectorial fuzzy variable (N > 1) described by a (joint) membership function μz. The support of μz is defined as For α ∈ [0, 1], the α-cut C α of μz is defined as If the joint membership function associated withz has the form μz = min i μz i thenz is called non-interactive and interactive otherwise.
Fuzzy sets are quite versatile since there are no restrictions on the (type and structure of the) used sets. However, for numerical methods to become efficient, certain assumptions are beneficial. The first restriction mimics numbers and vectors.
Definition 2 (Fuzzy number/vector) Letz ∈ F(Z ) with Z ⊂ R N for some N ∈ N such that Z is bounded and convex and the (joint) membership function μz is upper semi-continuous, i.e. lim sup and quasi-concave, i.e.
If there exists a unique z * ∈ Z such that μz(z * ) = 1 then we callz a fuzzy number for n = 1 and a fuzzy vector for n > 1.
This notion is easily extended to intervals and domains.

Definition 3 (Fuzzy interval/domain)
With the same assumptions as in Definition 2, there exists a subset S ⊂ Z such that μz(z) = 1 for all z ∈ S. Thenz is called a fuzzy interval for n = 1 and fuzzy domain for n > 1 .
Note that the quasi-concavity of the membership function implies convexity of any α-cut C α . In particular, given S ⊂ Z as in Definition 3, the convex hull conv(S) is a proper subset of C 1 , too. The quasi-concavity also leads to the nestedness property of α-cuts from Definitions 2 and 3, i.e.
This property is essential for the α-cut propagation method, see Theorem 2. In the following example, we introduce the most common fuzzy structures. Namely, the fuzzy trapezoidal interval which we will use throughout this work.
The triangle fuzzy numberz = , z * , r specified by left and right limit , r and peak position z * = z * = z * r is a special case of the trapezoidal fuzzy set. Figure 2 depicts the propagation of a trapezoidal fuzzy set.
For arbitrary fuzzy sets Zadeh's extension principle is the way to propagate uncertainty through a mapping.
If more underlying structure is given, the extension principle can be formulated equivalently in terms of constrained optimization. Zadeh's principle can be reformulated into the so-called α-cut propagation for fuzzy vectors and fuzzy domains to reduce the computational costs. Theorem 2 (α-cut propagation [30]) Let f : Z → V be continuous between metric spaces (Z , d 1 ) and (V , d 2 ) and letz ∈ F(Z ) with C 0 [z] ⊂ K ⊂ Z for a compact set K with convex Z . Furthermore, let the membership function μz be quasi-concave and upper semicontinuous. Then μf can be characterised via α-cuts as Note that (9) follows from the assumption that C α [z] is closed and consequently compact as it is a closed subset of a compact set K in a metric space. Since f is assumed to be continuous, compact sets are mapped to compact sets.
Comparing Zadeh's principle with the α-cut propagation, we immediately see that the former requires the inverse image and an optimization step for each point v whereas the later only depends on two optimization steps of f for the number of α-cuts. Since it is infeasible to perform Zadeh's principle for all points v ∈ V , it is combined or replaced with a sampling approach. We distinguish two variants to sample in V : • Semi sampling approach: directly solve the constrained optimization problem with a global optimiser. For a given sequence (v k ) k ⊂ V , compute the supremum over Z k := {z ∈ Z : f (z) = v k }, see the red line in Fig. 2.
Use the data sample pairs (v k , μ k ) and reconstruct μf , e.g. by convex hull or an envelope approach, see the orange/purple graphics in Fig. 2.
If f has high evaluation costs, the propagation inevitably becomes difficult and costly. This for example is the case if f represents a finite element solution with a large number of elements. Thus, to still render the propagation feasible the number of evaluations of f needs to be computationally feasible. In case of applicability of α-cut propagation, the number of evaluations depends on the (global) optimization method. This number is possibly lower than the sampling counterparts. However, if the requirements of Theorem 2 are not met one may fall back to the sampling methods. By introducing a surrogate model f h ≈ f , the computational costs can be further reduced for both strategies. In essence, this approach can be interpreted as a reduction of evaluations of f , determined by the number of evaluations needed to obtain f h , assuming that the evaluation costs of f h are negligible.

Upscaling of heterogeneous linear elastic material
In the following we develop a numerical homogenisation method which yields a macroscopic material based on microscopic properties. The goal of the approach is to dispose of the computationally involved microstructure and construct an upscaled material surrogate with similar homogenised behaviour on a larger scale. Equivalent terms for "homogenised behaviour" are macroscopic, effective or upscaled behaviour [31], where the examined "behaviour" is subject to some quantity of interest (e.g. average displacement or stress) based on the system response.
In classical (asymptotic) homogenisation, an effective property is calculated based on the assumption of an infinite periodic domain [32]. The local microscopic structure is defined in terms of a representative volume element, which in our setting would consist of a single inclusion that has identical shape in the entire domain. For problems in nonperiodic media, the methods of numerical homogenisation or numerical upscaling as e.g. described in [33] are used, where local boundary value problems are solved to calculate effective characteristics in each local domain [34], also see [9,12,35] for analytical stochastic homogenisation. Periodic boundary conditions (PBCs) are commonly used for numerical upscaling methods of matrix composite material [36,37]. In this work, we focus on numerical homogenisation with PBCs for linear elasticity. Let D = [−1, 1] 2 be a unit cell domain and let the heterogeneous material law be encoded such that σ · n is antiperiodic on D with n denoting the outer normal with respect to ∂ D. Let and σ denote the average strain and stress of the computed strain and stress σ . Then, there exists a tensor H ∈ R 2,2,2,2 satisfying σ = H : .
The tensor H is called the effective or upscaled (macroscopic) tensor, representing the elastic moduli of the homogenised medium. By construction it holds that E = . Consequently, H can be obtained from (12) by choosing macro strains E corresponding to different elementary load cases and subsequently computing σ by solving (11). In what follows, the homogenisation technique will be applied for every shape parametrisation encoded in C.

Measuring the distance between constitutive tensors
holds. A tensor C ∈ Ela(d) may exhibit further symmetry properties. If there exist Lamé constants λ, μ ∈ R such that the tensor is called isotropic. For orthotropic tensors, we define a rotation matrix and the full orthogonal group operation o(g) via (o(g)C) i jk = g i p g jq g kr g s C pqrs , C ∈ Ela(d).
We then say that C ortho ∈ Ela(d) is orthotropic if it can be represented as for some isotropic and forth order tensors O 1 and O 2 with a Kelvin-Mandel matrix representation [38] given as With that representation, any orthotropic tensor has 5 degrees of freedom (λ, μ, 1 , 2 , θ). The orthotropic tensors represented in normal form have only 4 degrees of freedom since θ is chosen as zero. Let Ela(d, iso), Ela(d, ortho) ⊂ Ela(d) be the set of all isotropic and orthotropic tensors, respectively. Given any anisotropic tensor C, we define the distance to the symmetry class of isotropic and orthotropic tensors by with the Frobenius norm · . The distances to the symmetry classes can be characterised as follows.
Proposition 3 (Distance to isotropy class [39]) as the orthogonal projection of C onto the isotropic material class. Then, Proposition 4 (Distance to orthotropy class [39]) Then, With the introduced material class concepts of constitutive tensors we are able to measure the distance of our effective material to the isotropic respectively orthotropic material class for each shape parametrisation.

Fuzzy edge detection
Given a blurred image as depicted in Fig. 3, the most common approach to reconstruct the original image is to use some method from the wide class of blind deconvolution methods [40]. These methods assume that an image y is the convolution of an original image x and some kernel k, distorted with additive noise n, i.e., Finding a pair (x, k) satisfying equation (17) is equivalent to de-blurring the image, leading to a original image x. This however is an ill-posed problem, since an infinite number of such pairs can be found [41], making some form of regularisation inevitable. The classical approaches assume zero noise and employ regularised least squared methods [42][43][44]. Each regularisation is based on assumptions of the kernel and the original image. More recently, natural image statistics and the Bayesian framework were used to formulate such assumptions more precisely and thus improved the performance of blind deconvolution methods, see [45][46][47] for a first overview. Simply put, prior knowledge and the ability to incorporate this knowledge improve the deconvolution result. However, we are in a different situation since we miss this prior knowledge and we are shall not be inclined to make These three assumptions implicitly define a possible set of boundary curves provided a blurred image. Any of these boundary curves is possible but nothing can be said about the probability of each curve. Since we assume a star-shaped inclusion, a radial function r (θ ) : [0, 2π) → R + is able to describe the boundary of this inclusion. The uncertainty of the boundary is then captured in a fuzzy functionr .
In order to construct a fuzzy functionr from a blurred image with domain D = [−1, 1] 2 , we use trigonometric interpolation on a vector of N ∈ N fuzzy numbers. Each component of the fuzzy vector denoted bỹ is constructed from a radial cut (i.e. a line segment) that starts in the center and ends at the outer boundary for a prescribed set of angles determines the intensity of the gradient along this radial cut, see Fig. 4 for an illustration. The information encoded in I i is used to construct a trape-  (6) obtained from a intensity function I i for given angle θ i and thresholds 0 < δ < δ < 1 zoidal fuzzy interval. For this, consider lower and upper percentage threshold values 0 < δ < δ < 1 and determine M := max r ∈[0,R max i ] I i (r ). If the intensity in one point is larger than δ M, it is assigned a possibility of one. If it is smaller than δ M, it gets assigned a possibility of zero. For simplicity, in between these thresholds we assign a possibility by linear interpolation. This algorithm yields a membership function μr θ i of trapezoidal form, which is illustrated in Fig. 4.
Gathering these fuzzy numbers into a vector yields the non-interactive fuzzy vectorR N , encoded in the membership function where non-interactivity means that the value of one sample does not influence the membership function of another value. Let R = (r 1 · · · r N ) ∈ C 0 [R N ] and denote the trigonometric interpolation by It represents the mapping from interpolation points onto the trigonometric polynomials. This mapping is bijective if the degree of freedom coincides with the number of interpolation points. The coefficients are efficiently computed via a discrete Fourier transformation and the interpolation scheme yields convergence rates as follows.
Given all possible interpolation points R ∈ C 0 [R N ], we define the set of radial boundary functions describing the interface as Note that T (·;R N ) defines a fuzzy function [29]. The number of cuts and thus the number of interpolation points N is setted according to the accuracy required by the considered problem. With sufficiently many interpolation points, the set B N can represent highly oscillatory boundaries. For the sake of efficiency, the interpolation should be carried out with as few points as possible, which demands an adequate knowledge about the properties (in particular smoothness/roughness) of the physical system at hand. We choose the total variation of the radial function as a measure for roughness, which is Given some bound 0 ≤ b ≤ ∞, the TV restricted set of radial boundary function based on N grid points is defined by With this construction, we may define the interactive fuzzy set where 1 X denotes the characteristic function on a set X . In Sect. 3.1 we make clear that the fuzzy set in (21) in fact defines a fuzzy vector. This in turn motivates the αcut propagation of fuzzy uncertainty, introduced in Sect. 2.1.
In particular, consider a continuous real valued function Q defined on C 0 [R N ] which is the 0-cut of the non-interactive fuzzy vectorR N from (18). Note that C 0 [R N ] is compact since the image domain D = [−1, 1] 2 is bounded. The propagated uncertainty of the interactive fuzzy setR N ,b through Q is then captured on each α-level of the fuzzy setQ b := Q (R N ,b ). We emphasise that C 0 [R N ] defines a tensor domain on which Q is well defined, i.e.
a fact that becomes useful when applying low-rank tensor formats in the upcoming Sect. 4.2 as surrogate models for Q. However, the computation ofQ b only requires evaluation of , which in general is not a tensor domain.

Properties of the TV bounded fuzzy setR N,b
The total variation determines the shape of the inclusion.
With a very small total variation, the shapes become more and more circular. Any circular shape, independent of the radius, has a total variation of zero. We want to point out, that we do not measure the total variation of the trajectory, where the total variation of a circle would be larger than zero. Instead we measure the total variation of the radial function. In Fig. 5 shapes with different total variations are depicted. The blue shape is generated by taking alternating radii, the white shape is the result of optimization. The maximal total variation is 5.2 for N = 12, whereas a randomly generated shape has a total variation of around 2. With increasing N the maximal total variation would also increase.
Note that the resulting shapes may violate the boundaries in between the sampling points. Especially, the maximum total variation solution. This can be resolved with more sampling points, a sufficiently strict TV bound or by replacing the trigonometrical interpolation with a corresponding spline interpolation.
Bounding the total variation leads to interaction of the fuzzy setR N ,b . Figure 6 illustrates the interaction for the case N = 3. It shows that fixing one point constraints the remaining points to the light gray area. Hence the possibility outside this area is zero. In a non-interactive setting the possibility could be strictly larger than zero. Consequently, the total variation bound restricts the set of possible curves. The Proof Let R 1 , R 2 ∈ suppR N , t ∈ [0, 1] and define R t = t R 1 +(1−t)R 2 . The radial function to describe the boundary then takes the form Consequently the α-cuts ofR N ,b are nested and convex. Note that for given N , there always exists b = b(N ) such that the conditions of Proposition 7 hold true. In particular, for b large enough it holdsR N =R N ,b . Consequently, the α-cut propagation (23) can be applied.

Remark 1
The convexity property from Proposition 6 is essential to apply Theorem 2. It is a consequence of the linear interpolation and the triangle inequality property of the TV based restriction. In the case, one of the latter properties does not hold, the fuzzy set may not define a fuzzy domain. Then, the propagation of the fuzzy set still can be realised by sampling approaches as discussed in the end of Sect. 2.1.

Accelerated emulation of fuzzy effective material
The fuzzy edge detection described above yields a fuzzy vec-torR N , where N is the number of radial cuts. Trigonometric interpolation for each realization R ∈R N then results in a boundary representation of the inclusion. In the following, we define a composite material based on this representation. Consider the square [−1, 1] 2 on which we define a two phase composite material C(x) = C(x, R) with C incl , C matrix ∈ Ela(d, iso). This description represents a piecewise isotropic material with star-shaped inclusion defined as Given a fixed R, the homogenisation from Sect. 2.2 yields the constitutive tensor C(·, R) ∈ Ela(d) for the upscaled macroscopic material. In Voigt notation, this mapping is denoted as with symmetric positive definite matrix H(R) ∈ R 3,3 . It in general describes an anisotropic material since the involved geometry of the inclusion may lack any type of symmetry. We would like to point out that the parametric dependency R → C(·, R) does not define a continuous function with images in L ∞ (D) d,d,d,d since marginal changes of the shape of the inclusion immediately yield an L ∞ error equal to the contrast C incl − C matrix F with Frobenius norm · F . Despite this discontinuity and its effect on the regularity of the parametric solution R → u(R), the parametric homogenised tensor (26) defines a continuous map.
Recall that the evaluation of H involves multiple simulations of periodic linear elastic boundary value problems of the form (11). To accelerate the upscaling process, in the following we replace the simulator with an emulator. The emulator relies on a mesh discretisation of the domain D and a group of tensor train surrogates to approximate H, which is discussed in Sect. 4

.2.
If D is discretised by an automatic mesh generator under the constrained of equal amount of vertices on the inclusion's boundary and on the boundary of D for any R, the resulting mesh map R → M(R) in general is discontinuous. As a consequence, a mesh dependent finite element computation may inherit the lack of continuity even through the whole geometry has smooth dependence on R, see Fig. 7 for an example. This in turn aggravates optimization based on gradient information. We solve this issue in Sect. 4.1 by constructing a family of transformed meshes with smooth dependence on R.

Constructive smooth transformation of meshes
For n ∈ N interpolation points, consider the composite interface discretisation We now construct a smooth diffeomorphism that creates meshes based on a single reference meshM denoted as Compute the sets F 2n (R ref ) and F 2n (R) corresponding to higher discrete resolution of the reference or target interface. The transformation then consists of two parts as follows. First, let {φ i } form a smooth partition of unity of [0, 2π ] with φ i (iπ/n) = 1 and supp φ i = [(i − 1)π/n, (i + 1)π/n]. This function set deals with the angular part of the transformation. Second, let χ i ∈ C k [0, √ 2] be a set of splines with k ≥ 2 for the radial part of the transformation with the following properties: 1. The reference radius is mapped onto the transformed radius, i.e.
2. The spline is strong monotonically increasing on (R min , R max ) and otherwise equals the identity map.
Altogether, the transformation reads Note that the number of radial cuts N determines the range of shapes, whereas the number of interpolation points n has to be sufficient large to trace the shape and to account for the resolution of the finite element mesh. We refer to Fig. 8 for an illustration of the capacity of the constructed map , which is characterised next.  Proof The bijectivity follows from strong monotonicity in radial direction and the smooth partition of unity in angular direction. The smoothness of φ i and of the polar coordinate mapping away from zero is given since χ i is k-times continuous differentiable. For i = 1, . . . , 2n the desired property follows immediately.

Tensor trains based emulation
Assume some function f : based on a polynomial feature class . . , N and unknown coefficient tensor U : → R. This model suffers from the curse of dimensionality since the cardinality | | grows exponential with respect to the dimension N . A possibility to circumvent this challenge lies in a compressed representation of the coefficient tensor, here based on the tensor train format (TT format) [48]. Let ρ = (ρ 1 , . . . , ρ N −1 ) ∈ N N −1 be the tensor train rank and let ρ 0 = ρ N = 1. We then choose U given in tensor train decomposition by with component tensors of order 3, such that U i [:] ∈ R ρ i−1 ,| i |,ρ i for i = 1, . . . , N . The number of degrees of freedom for this design is given as N i=1 ρ i−1 | i |ρ i , which does not grow exponentially but only linearly in dimension N . Given data (R (k) , y (k) ) K k=1 , K ∈ N with y (k) = f (R (k) ), we obtain a surrogate as in (28) by carrying out a regularised empirical regression, namely by solving the optimisation problem min U as in (29) K k=1 f (R (k) ) − y (k) + λ U Fro (30) with regularisation parameter λ > 0 and the Frobenius norm · Fro . The optimisation problem (30) can be solved by regression with an alternating linear scheme (ALS) [49]. As a modification of the basic ALS, we introduce a scheme for rank adaptivity. This concept offers several advantages to obtain an adjusted tensor train rank which is in general unknown a priori. From a practical point of view, it reduces the computational cost during optimisation while obtaining a prescribed accuracy in the approximation class. Furthermore, starting with a small rank, the iterative process empirically enables the ALS to converge to a solution based on successive computations of initial guesses for models with higher ranks based on the given restricted number of samples. The proposed approach is summarised in Algorithm 1 and the rank adaptivity is presented in Algorithm 2.
A principle of the rank adaptivity is to keep one (control) singular value to monitor the importance of the related rank coupling value during the optimisation process. This additional singular value per rank coupling remains until the end of the regression scheme to prevent oscillation between rank growth and reduction. It ensures an upper bound of the related rank value throughout the entire process. After successful termination of Algorithm 2, the existing (possibly small) control singular value is removed by a final rounding [48] of the resulting tensor train. We refer to [49] for more technical details on the basic ALS, e.g. with regard to setting and moving the non-orthogonal component (the core) via left and right matricisation and orthogonalisation. Eventually, the overall emulator is obtained by evaluation of the 6 scalar valued tensor train surrogate maps, replacing R → [H(R)] i j and 1 ≤ i ≤ j ≤ 3 corresponding to the upper triangular part of the symmetric matrix H(R).
Note that non-zero values of H 13 and H 23 solely appear when the effective material behaviour is not isotropic or rotational-free orthotropic, i.e., this is described by (13) with θ = 0.

Algorithm 1 Rank adaptive empirical tensor train regression via ALS
Require: Define operators L and R as left and right matricisation of order 3 tensors.

Numerical experiments
This section is concerned with the assessment of the numerical performance of the approach presented above, for which two experiments are examined. In the first experiment, we measure the distance of the upscaled material to the isotropy and orthotropy class. Additionally, we identify configurations that maximise and minimise the respective distances.
In the second experiment, we test if the upscaled material is suitable for a worst case analysis, i.e., we test if it is possible to find bounds that envelop the most extreme behaviour of the quantity of interest under consideration. For this, a material with 8 × 8 arbitrary inclusions placed on a checkerboard is compared to an upscaled material with the same geometrical dimensions. Both experiments demonstrate the efficacy of the fuzzy approach to model uncertainty of the inclusion to identify extremal behaviour and its source. Three types of computational tasks are performed in the experiments. The foundation is laid by Finite Element (FE) simulations to generate the realisations of constitutive tensors by solving (11). Furthermore, an alternating linear scheme is used to train the surrogate model and for optimization computations to carry out the uncertainty propagation based on α-cuts.
All FE computations are done with the FEniCS package [50]. The mesh generation is realised using Gmsh [51] and the Bubbles package [52]. Moreover, the python package TensorTrain [53] is utilised for the rank adaptive tensor train regression and ALEA [54] for the underlying polynomial features. The optimization tasks are performed with restarted trust region optimization implemented in Scipy [55].
The computations are performed with N = 6 and N = 12 radial sampling points. Here, the N = 6 case imposes significantly less computational burden on the surrogate and the optimization than the N = 12 case. In particular, our optimization scheme took up to 10 6 evaluations of the H(R) for N = 12 and up to 1 × 10 5 evaluations for N = 6 per propagated quantity of interest. We underline that all optimization tasks involved are non-convex with non-linear cost functions and constraints. Consequently, we take arbitrary points R ∈ [0.3, 0.7] 6 for N = 6 and R ∈ [0.4, 0.6] 12a smaller domain-for N = 12. This enables sufficiently non-trivial shapes without rendering the surrogate model and the optimization infeasible. Each realisation R determines a boundary stellar inclusion in terms of a trigonometric interpolation. Inside and outside of the inclusion we assume isotropic material behaviour. Concretely, we let the Young's modulus E incl = 3230 and Poisson ratio of ν incl = 0.3 in the inside. Outside of the inclusion we choose E matrix = E incl /4 and ν matrix = 0.2. These values are transformed into Lamé constants by By adaptation to the plane, adapted Lamé constants are obtained, namely For a given R ∈ C 0 [R N ] = [0.3, 0.7] N the homogenisation problem (11) is solved with FE of uniform polynomial degree p = 2. The reference mesh is based on 80 nodes on the related reference composite interface T (·; R ref ) is shown in Fig. 8. The computational mesh with composite interface T (·, R) is then obtained by applying the transformation (·, R).

Surrogate validation
As a preparative step, we build an accelerated emulator in the tensor train format for N = 6 and N = 12 dimensional fuzzy input vectorsR N , yielding a compressed surrogate of the map R → H(R).
Given K = 15635 (N = 6) and K = 12000 + 4096 (N = 12) samples that are normalised with respect to the sample mean and variance, we iteratively apply Algorithm 1 and Algorithm 2 with tensorised Chebyshev polynomials up to degree 5 in each coordinate. More precisely, we train the initial tensor train surrogate with a degree of two. Next, a new tensor train for polynomial degree deg + 1 up to deg = 4 is set to the last tensor approximation as initial guess instead of starting with a random rank-1 tensor train. All tensor train component entries associated to the higher polynomial degree are initially set to zero. This training approach turn out to have two advantages. First, in the iterative procedure the training and validation sets are split randomly, resulting in a pseudo stochastic solver. A strategy that is successful in the context of machine learning. Second, there is a huge speed-up in the training procedure for finding a local minimum. In fact, with a naive tensor train training for N = 12 with initial tensor polynomial degree 5, we rarely observed convergence to meaningful surrogates given random initial values.
We use tol MSE = 1 × 10 −8 , iter max = 200, L HIST = 10, tol DECAY = 1 × 10 −5 , r max ≤ 5, 7, 9, 20 for deg = 2, . . . , 5 , r add ∈ {1, 2}, and δ = 1 × 10 −8 . The results of the surrogate learning approach are depicted in Tables 1 and 2, which complement Figs. 10 and 11, showing mean square error or pointwise relative and absolute errors. The highlighted entries in both tables show a pointwise lack of accuracy of the tensor train surrogate for the H 13 and H 23 components, which represent the pure rotational or anisotropic contribution in the overall upscaled tensor H on the used test sets. A further inspection of the high relative errors showed that these strictly belong to pointwise approximations of values close to zero. Nevertheless, strikingly small mean squared errors can be observed as a result of the proposed optimization algorithm. Figures 10 and 11 display the value ranges of the subcomponents of H and the overall absolute and rel-   ative errors of the subcomponents and the full tensor H. We observe that the reduced (pointwise) accuracy of the surrogates for H 13 and H 23 does not influence the overall error for the full parametric tensor H. The final obtained ranks obtained with the rank adaptive algorithm are depicted in Fig. 9 for N = 6 and N = 12.

Distances to the isotropy and orthotropy material class
In this section we apply the tensor train emulator to obtain an approximation of the parametric constitutive tensor H(R) of the effective material defined in (26), given an parametrisation of the inclusion through the vector For each constitutive tensor we compute the distance to the isotropy class, denoted by d iso , and orthotropy class, denoted by d ortho , according to (14). The fuzzy uncertainty of the boundary is propagated in terms of α-cuts. On each α-level, two optimization problems with constraints defined in (23) have to be solved, where Q is either d iso or d ortho . As opti-miser we use a restarted trust region scheme to minimise the non-convex and non-linear target function with non-linear constraints. Figure 12 shows the experimental results for Q = d iso with N = 6, zero α-cut C 0 [R N ] = [0.3, 0.7] N and the total variation bounds b = 0.5, 1, ∞. We point out that since D = [−1, 1] 2 is not of special hexagonal shape but of square shape, a circular inclusion does not yield isotropic effective material properties. Nevertheless, the circular shape-while exhibiting vanishing total variationprovides effective properties closest to the isotropy class. Since in our experiments we choose a stiffer material as inclusion, the circular inclusion minimises its size proportional to the possible maximal radius encoded in the α-level. In this sense, the size of a circle is a measure of perturbation of the isotropic matrix material. Moreover, if the total variation and radii bounds allow it, the inclusion converges to a peanut shape. The green peanut shaped inclusion in Fig. 12 marks the inclusion with maximum distance to the isotropic material class for the given parameterisation family B N ,∞ from (20). Note that by rotational invariance also the peanut rotation of π/2 yields the same maximising shape. Figure 14 pictures the result for Q = d iso with N = 12, zero α-cut C 0 [R N ] = [0.4, 0.6] N and total variation bounds b = 0.5, 1, ∞. The shrinked domain has an immediate impact on the maximal possible total variation bound of all interface realisations. Consequently, maximum distances are caused by curves with similar shape up to rotation. Already for a total variation of 0.5, a symmetric shape emerges that becomes more prominent with increasing total variation. Figure 13 depicts the experimental results for Q = d ortho with N = 6, zero α-cut C 0 [R N ] = [0.3, 0.7] N and total variation bounds b = 0.5, 1, ∞. A minimum of 0 is attained on each α cut of the corresponding membership function, since the orthotropy class Ela(2, ortho) contains the square symmetric class. Many boundaries in B N ,b lead to a shape that yields a square symmetric effective material.
As typical maximiser of the distance to the orthotropy class, a bean shape is found. We observe that the maximal possible total variation value of B N ,b is not exhausted when attaining the maximal distance. Note that due to rotational invariance, the bean shape can be rotated by multiples of π/2 while still being a (in fact the same) maximiser. Figure 15 shows the experimental results for Q = d ortho with N = 12, zero α-cut C 0 [R N ] = [0.4, 0.6] N and total variation bounds b = 0.5, 1, ∞. Again the shrinked domain has an immediate impact on the maximal possible total variation bound of all interface realisations. Consequently, maximum distances are attained by curves with similar shapes up to rotation. We observe that symmetric dented shape maximises the distance to the orthotropy class for the α = 1 cut or a total variation smaller or equal to 0.5. However, as the total variation bound gets larger and we allow for larger fluctuations due to a wider parameter range the inclusion converges to a non-axis aligned structure with non-symmetric dented shape.

Best/worst case estimates for non-homogenised matrix composites
Let On each subdomain D i j we consider a single two phase composite as defined by (25) using a local polar coordinate system w.r.t. the midpoint of D i j , encoded in a material tensor C inhomo with piecewise values C incl and C matrix specified in Sect. 5.1.  For this experiment we choose the non-interactive fuzzy setR N ,∞ with N = 6 trapezoidal fuzzy components characterised by In between these prescribed points, the α-levels are obtained by linear interpolation. The TV bound is ignored in this experiment. Furthermore, we consider a homogenised material tensor C homo = H(R) for R ∈ C 0 [R N ,b ] obtained from the homogenisation process of Sects. 2.2 and 5.1. In each domain, we model individual star shaped material inclusions as in (25) with boundaries of uniformly bounded total variation smaller or equal to b.
representing a linear elastic tensile test.
We are interested in the fuzzy propagation of the homogenised material through the average displacement Q = u defined by with the solution u of (31). Figure 16 shows the resulting membership function based on the fuzzy homogenised material law. The black dots mark resulting average displacements obtained for various non-periodic inclusions. Moreover, some configurations of possible realisations of Right: underlying coloured shapes of the inclusions associated with the extremal points displayed as square markers with same colour along the membership function. (Color figure online) the 8 × 8 inclusions are displayed. The corresponding average displacement is marked with colours, accordingly. Given any non-periodic configuration of prescribed matrix composites, we observe that the membership function at each α-level bounds the various average displacements. In particular, for α = 0., 0.5, 1 we illustrated a corresponding single composite shape that (assumed in a periodic structure) results in the minimum and maximum average displacements. As the interior composite is slightly stiffer than the surrounding matrix material, it is to be expected that the shapes of the minimiser and maximiser attempt to exhaust or avoid the maximum capacity of the parameter range encoded in C α [R N ,∞ ]. It is noteworthy that the homogenised material is able to serve as worst/best case estimator for this particular quantity of interest (i.e. the average displacement). As long as the repetition of the same cell leads to extremal values of the quantity of interest, the fuzzy homogenised material functions as a worst/best case estimator.

Conclusion
We considered the possibilistic shape uncertainty of one inclusion induced by a blurred image. In this setting, com-putational stochastic homogenisation was carried out to propagate the uncertainty through the linear elasticity model, resulting in a fuzzy effective material. Eventually, the effective material was examined with respect to its distance to the orthotropic/isotropic material class. Moreover, it was used as a worst/best case estimator with respect to the global average displacement for a non-homogenised material.
To achieve this result, two major obstacles had to be overcome. The first was caused by the discontinuous mapping from boundary to effective material, which was a consequence of the automatic re-meshing for each new boundary instance. This was successfully resolved with an arbitrarily smooth transformation of a single reference mesh onto the respective mesh for a domain with a star-shaped inclusion. The second challenge was given by the computational effort needed to perform the (possibilistic) uncertainty propagation. Depending on the optimiser and the quantity of interest at hand, each α-cut required up to 5 × 10 4 computational homogenisation simulations. This problem was successfully resolved with a highly accurate tensor train emulator.
We provide an (numerical) analysis for a problem where only little knowledge is given and only simple assumptions can be made. This may be distinguished from other works, where more prior knowledge is given and more elaborated constellations in C inhomo . Right: various periodic (red and green) and non-periodic (blue and magenta) composite constellations yielding average displacements within the area of black dots assumptions are suitable. In the latter scenario, one usually can utilise the Bayesian framework to model the uncertainty, see [56,57], where in our scenario stochastic modelling may lead to a false confidence [23]. While usually the uncertainty of material properties is modelled with precise probabilities, leading to statements about expectation values and probabilities of failure, cf. [58,59], we model geometrical uncertainties with imprecise probability, leading to statements about possible intervals and worst/best case configurations. Even if an imprecise probability approach for geometrical uncertainty is chosen, the geometrical shape is often restricted to simple parametrisations like a circle of varying diameter [60]. This in fact is a very useful simplification, for instance if the observed quantities only depend on the volume ratio of the two materials. Nonetheless, our experiments demonstrate that the shape determines if the effective material is isotropic, orthotropic or anisotropic. A circle always results in an orthotropic effective material in the considered case with unit square cells. Furthermore, our tensor train surrogate shows a remarkable accuracy of a relative error of order O (10 −4 ). This is partly the result of a suitable choice of features, the novel rank-adaptive training strategy and the continuity induced by using a mesh transformation instead of automatic re-meshing. Compared to other works such as [61,62], where standard generic methods like artificial neural networks and polynomial chaos are used in a similar context, we gain several orders of accuracy.
To focus on the influence of the shape uncertainty, other parts in our setting and analysis were kept simple. These parts are ideal starting points for future research. The homogenisation, for example, was performed with two phases of linear elastic materials. This can easily be substituted with more sophisticated materials-one may think of anisotropy, damage and higher contrast-or different homogenisation settings. The assumption of a star-shaped inclusion can be discarded with an adjusted edge detection and a different parametrisation of the boundary. Furthermore, the roughness measurement via total variation can be exchanged or extended with other restrictions of the boundary curves. In addition, it is possible to model the material constants appropriately with precise probabilities. By doing so, we add a stochastic dimension which together with the possibilistic uncertainty results in a fuzzy-stochastic model.
Although we only investigated 2D materials, our approach can be extended to 3D. The numerical upscaling and the surrogate model are applicable in higher dimensions. For the fuzzy edge detection the radial cuts should be replaced by cuts through the sphere. Under the assumption of star-shaped inclusions in 3D, the inclusion surface can be described in spherical coordinates with suitable interpolation for the radial component, which in turn can be constrained by a corresponding 2D TV measure. Moreover, in the 3D case there are more material classes. Depending on the symmetry class, the associated distances are-to our knowledge-not yet ana-lytically given and thus have to be computed numerically [39,63]. Thus, the theoretical adaption is feasible, but we expect the computational challenges of the 3D case to be significantly larger.
Another interesting research direction for our framework are more involved material models like hyperelastic materials or materials with damage modeling. In this case the numerical upscaling leads to an non-linear effective material law. In order to represent this non-linearity in our framework the surrogate model design needs to be adapted, e.g. mapping to interpolants of effective potentials based on [64].
Finally, we hope that our presented contributions, namely the treatment of discontinuity from re-meshing, the fuzzy edge detection, and the highly accurate tensor train surrogate, can be used and extended beneficially in other research to simulate composite materials in the presence of shape uncertainty.