Component Mapping Automation for Parametric Component Reduced Basis Techniques (RB-Component)

The aim of this paper is to develop some techniques for automation of the mappings (between working and reference domains) required by reduced basis methods: the development of geometry mappings is indeed often a substantial impediment to the implementation of reduced basis techniques, especially in the context of the reduced basis element method (RBEM) and the reduced basis component method (RBCM). In the RBCM context, the geometry mappings are applied at the level of components. The methods have been tested on various cases to understand the limits of the approach and try to foresee and overcome the possible failures


Introduction and motivation
The objective of this project is to compute -in real-time -the solution to parameter dependent partial differential equations (PDEs), where the parameters include geometrical factors generically denoted here as ρ.In this paper, the PDE we consider as an example is the Laplace problems (1) set on the spatial domain Ω ρ ⊂ R d (d = 2 or 3) with varying values of given sets of geometrical factors ρ : Find φ ∈ H 1 (Ω ρ ) such that where the boundary ∂Ω ρ is composed of two parts: ∂Ω ρ = Γ f ∪ Γ ρ with Γ ρ the parameter dependent boundary of the parametrized domain and Γ f the fixed boundary (possibly empty).
[-] and g f is an appropriate functions.Classical discretization techniques, such as finite element methods, may be too expensive if multiple resolutions are required or real-time response is expected.In this perspective, the Reduced Basis (RB) method [1,16,18] exploits the parametric structure of the PDE to construct fast and computationally efficient approximations.
To be even faster the RB method may be combined with Domain Decomposition and leads to componentbased RB approaches namely the reduced basis element method (RBEM) [14,15] and the reduced basis component method (RBCM) [11,12,19].In these versions of the reduced basis method, the domain of interest Ω ρ (where the PDE is set) is decomposed into a series of subdomains with simple shapes called components Ω ρ = ∪ K k=1 C k,ρ [5,10].Let us consider for example the case of Fig. 1 that represents the spatial domain of a horn for which the length L and the radii a 0 and a mouth can vary [13].
Figure 1: The decomposition that is proposed in the frame of the RBCM is exemplified on Fig. 2. Each of the components featured in there is obtained by deformation of one reference component chosen among a set of few reference components.Each reference component is provided with some basis functions (reduced basis functions) that represent the behavior of the set of all the PDE solutions on such subdomains.The restriction of the solution to (1) to every component is then sought as a linear combination of those basis functions mapped onto the component from the associated reference component.
The objective of this RB-Component project is to rapidly propose the mapping that needs to be used to transfer back and forth all the informations (mesh, reduced basis, geometrical factors) from the reference components to each associated subdomain in Ω ρ in order to solve the PDE of interest on the global domain Ω ρ .
It is not uncommon to use the elasticity equation to lift boundary data into the interior for the purpose of geometry mappings (for example, for ALE fluids calculations) or mesh generation, see [21,20,8,6].
In our approach we propose new strategies to generate these maps using the solution u u u to a linear elasticity problem.Each new subdomain is obtained by deforming the appropriate reference component through a map T : x x x → x x x + u u u(x x x) that only has to send the boundary of the reference component onto the boundary of the subdomain.
The different methods that have been tested differ from the way the displacements u u u are imposed on the component boundary.The first technique that we have investigated, consists in simply imposing classical Dirichlet boundaries conditions, but it requires an explicit parametric definition of the boundary component, which is not always possible.The second approach that we have studied, consists in a penalization method which only requires an implicit caracterization of the boundary; though this description of the boundary is less precise it appears sufficient for our purpose.
In what follows we shall focus on a single component, in the sense that we do not perform any domain decomposition.There is thus only one reference component and the corresponding deformed component.In Sec. 1, we tested these different approaches to build an automatized mapping in order to solve Laplace problems as (1) over different kinds of geometry.Then, in Sec. 2 we present the mapping used in a reduced basis method context.) the coordinates of the associated point x x x in Ω ρ .We choose a mapping T T T as follows where u u u(x x x) = (u 1 (x x x), u 2 (x x x)) is the displacement.As said above, we choose that the displacement is the solution of a linear elasticity problem over the reference domain 2µ e e e(u u u) : e e e(v v v) where V is a space that is in some sense defined as (the precise definition being given hereafter) with the notation e e e(u u u) : e e e(v v v) := i,j e e e i,j (u u u)e e e i,j (v v v).
It is classical, and will be instrumental for our approach, to remind that problem (3) is linked to the problem of minimizing the energy e e e i,j (w w w) 2 )dx x x. ( and the first choice of space V , in line with this minimization process, is called here an explicit version (called also "pointwise") : assuming that s → x x x(s) (resp.s → x x x(s)) is a one-to-one parametrization of Γ ref (resp. of Γ ρ ): At this level it is interesting to recall that there are two ways to impose Dirichlet boundary conditions: the strong one where the discrete solution belongs to V and the weak one where the boundary condition is satisfied through a penalization formulation.In what follows we will first test these two ways and then focus on the weak one that appears much more simple to implement.In addition the weak formulation only requires an implicit caracterization of the boundary which is generally much more simple than having a parametrization (especially in higher dimension that will be dealt in a future paper Ω ρ ∈ R 3 ).This leads us to introduce an implicit version (called also "slippery") of the space V : assuming that Γ ρ is defined as the set of points x x x in R 2 such that F ρ (x x x) = 0: Let us now proceed to the use of the map from the reference domain: a simple change of variables leads to where J is the Jacobian matrix of T T T : and J is det(J ), more precisely, this reads One convenient way to verify that the mapping as defined above is correct, with respect to our aim which is to simulate partial differential equation on Ω ρ , is to consider the Laplace problem (1) over Ω ρ that we state here under a weak formulation: In what follows we tested different approaches to compute the displacements u u u which will be used in the mapping T T T for several test cases.All the simulations have been done using P 1 finite element within Freefem++ [7].Several examples are presented to compare different solution algorithms but also different treatments of F ρ (e.g., explicit versus implicit).More precisely, we first present in Sec.1.1.1,a circular hole with pointwise Dirichlet strong boundary conditions, which are the natural way to impose the displacement on the boundary Γ ρ .As it is difficult to use strong Dirichlet conditions for the general case, and we propose instead to use a penalized approach: we provide in Sec.1.1.2a comparison with penalized Dirichlet conditions.Then we consider a uniform shear square (with linear F ρ ) with strong (in Sec.1.1.3),and penalized (in Sec.1.2.1)Dirichlet conditions.The conclusion is that we do not loose much when using a penalized formulation.Thus we use a penalized version to treat increasingly more complex deformations: a non-uniform shear square (with nonlinear F ρ ) in Sec.1.2.2, then a bell (with nonlinear and large amplitude F ρ ) in Sec.1.2.3, and finally a crescent moon with a cusp, in Sec.1.2.4.

Example 1: square with a circular hole and strong Dirichlet boundary conditions
In this example, we consider a generic component in the form of a square [−1, 1] × [−1, 1] where a disc of radius a and centered in c = (0, 0) has been removed, the reference will have about the same shape except that the removed disc will have a radius equal to ā (see Fig. 3).In this example the set of parameter ρ is made up of the radius a.In the generic parameter dependent domain Ω ρ , the boundary Γ a represents the parameter dependent boundary Γ ρ .To compute a displacement u u u that describes the mapping from the reference domain Ω ref to the generic domain Ω ρ , we solved a linear elasticity problem with homogenous Dirichlet boundary condition on the fixed boundary Γ f , and non homogenous Dirichlet boundary condition on Γ ref -which represents the boundary of the reference domain that will be deformed in order to get the generic boundary Γ ρ = Γ a : In what follows we choose to set the Young modulus to E = 1 and the Poisson ratio to ν = 1 4 .
Figure 4: reference mesh T ref (left), deformed mesh T map (middle) and true mesh T h of the generic domain (right).
In Fig. 4 we represented the mesh T ref -a regular triangulation with 200 vertices on Γ a and 50 vertices on Γ ref -associated to the reference domain (for ā = 0.2) (left), and meshes associated to the generic domain (for a = 0.3): the deformation T map of T ref due to displacement u u u (middle) and the true mesh T h of the generic domain Ω ρ (right) built independently but similarly as T ref .We observe that the deformed mesh fits well with the objective (larger red circle represented on left and middle meshes of Fig. 4) and, at least from the eye point of view, appears rather regular and quite similar to the mesh T h with 200 vertices on Γ f and 50 vertices on Γ a .
In order to better quantify the quality of a mesh, we produce in Table 1 below classical quantities associated to the mesh.As usual, we denote by h, h min , h mean the maximum, minimum and average mesh size of T h , respectively.We also introduce σT = h T ρ T , where ρ T is the diameter of the incircle of a triangle T ⊂ T h , and for reference σ T = σT σ T , where T is an equilateral triangle.Then σ, σ min , σ mean denote the maximum, minimum and average of σ T for T ⊂ T h , respectively.These results highlight the regularity of the mesh obtained by our transformation.We denote by X 0 h the P 1 finite element approximation of H 1 0 (Ω ρ ) associated to the mesh T h and by X h the P 1 finite element approximation of X = {z ∈ H 1 (Ω ρ ); z = f on Γ f ; z = 0 on Γ ρ } associated to the mesh T h .Let φ h ∈ X h be the solution of the true discrete approximation of ( 8)

Mesh h
We denote by X 0 map the P 1 finite element approximation of H 1 0 (Ω ρ ) associated to the mesh T map and X map the P 1 finite element approximation of X associated to the mesh T map .Let φ map ∈ X map be the solution of the following discrete approximation of ( 8) where K is the mapping matrix defined by (6).
In order to validate this mapping approach to solve our Laplace problem (8) with g f = 1 and Γ ρ = Γ a .We have computed a true discrete approximation on the mesh T h and approximation using the mapping approach with ā = 0.2 and a = 0.3.
In Fig. 5 we show the solution φ map (left) and the relative error measured in the L ∞ -norm between φ h and I h φ map , where I h is the interpolation operator from X map into X h .This error is on the order of h 2 , where h is the maximum mesh size of T h (see Table 1), as might be expected from usual interpolant estimates, the order-unity derivatives for the data given, and the geometric factors.This boundary condition is imposed in a weak form, using a (quadratic) penalty method: we replace the constrained minimization problem inf w w w∈V0; w w w−(ā−a)n n n=0 by the unconstrained problem inf with and Let u u u be the solution of minimisation problem (12) and u u u the solution of minimisation problem (13), we have Besides, finding a solution to the minimisation problem ( 13) is equivalent to finding a solution to the following variational problem: which can be rewritten as follows : Find u u u ∈ V 0 , such that ∀v v v ∈ V 0 , 2µ e e e(u u u) : e e e(v v v) + λ div (u u u) div (v v v) dx x x + 1 Figure 6: Left: reference mesh T ref ; Right: deformed mesh T map .
In Fig. 6 we represented the mesh T ref associated to the reference domain (for ā = 0.2) (left), and the mesh associated to the generic domain (for a = 0.3) : the deformation of T ref due to displacement u u u (right), which is very similar to the mesh obtained in the previous approach where we imposed non-homogenous Dirichlet boundary conditions on Γ ρ = Γ a .As done in the previous example, in order to better quantify the quality of a mesh, we produce in Table 2 quantities associated to the mesh T map obtained by our mapping (see Sec. 1.1.1).The quantities associated to T ref and T h are already given in Table 1.Again, the results highlight the regularity of T map .Moreover, we observe that the results of   We consider the same Laplace problem as in the previous example.In Fig. 7 we show the solution φ map (left) and the relative error measured in the L ∞ -norm between φ h and I h φ map , where I h is the interpolation operator from X map into X h .This error is on the order of h 2 (see Table 2), as might be expected from usual interpolant estimates, the order-unity derivatives for the data given, and the geometric factors.

Example 3: deformed (uniform shear) square with strong Dirichlet boundary conditions
In this example we consider the deformation of the unit square [0, 1] × [0, 1] as in Fig. 8.

Γ3
ref The set of varying parameters ρ is made of the coefficient α and β.In the generic domain Ω ρ , the parameter dependent boundary Γ ρ is made of the union of Γ i ρ , 1 ≤ i ≤ 3 and respectively Γ ref is made of the union of Γi ref , 1 ≤ i ≤ 3.In order to compute the displacement that describe the mapping from the reference domain Ω ref to the generic domain Ω ρ we associate the following Dirichlet boundary conditions to the linear elasticity problem In Fig. 9 we represented the mesh T ref -a regular triangulation with 50 vertices on Γ f and each Γ i ρ , 1 ≤ i ≤ 3 -associated to the reference domain (left) and the deformed mesh T map for α = 2 and β = 1 (right).As previously, we observe that the mapped mesh fits well with the objective (in green).We consider now the Laplace problem (8) with g f = (1 − x 2 )x 2 .In Fig. 11 we show the solution φ map (left) and the relative error measured in the L ∞ -norm between φ h and I h φ map , where I h is the interpolation operator from X map into X h .This error is approximately O(h 2 ) where h is the mesh size of T h .As expected, the error is on the order of h 2 as in the previous examples.

Computation of the displacement using a penalty method
The technique that we have investigated in section 1.1 consisted in simply imposing Dirichlet boundary conditions in order to control the displacements u u u on the component boundary such that the deformation of the reference domain Ω ref matches with the generic domain Ω ρ .However, this approach requires an explicit parametric definition of the boundary component, which is not always possible.The second approach, that we now present, only requires an implicit caracterization of the boundary Γ ρ .This is done by the use of a functional F ρ defined such that F ρ (x x x) = 0, on Γ ρ .
The idea is to compute a displacement u u u such that F ρ (x x x + u u u(x x x)) = 0 on the boundary Γ ref , which leads to the following constrained minimization problem inf w w w∈V0, in which J(w w w) and V 0 are respectively given by ( 5) and ( 14).Nevertheless, we decided to weakly impose the constraint F ρ (x x x + u u u(x x x)) = 0 using a penalty approach, which leads to the following unconstrained problem inf w w w ∈V0 J(w w w ) + 1 which leads to solving the following variational problem: In what follows we consider different examples with affine or nonlinear function F ρ .In the nonlinear case, we propose different approaches to solve the problem, using a fixed-point method or a steepest descent method.

Example 4: deformed (uniform shear) square
In this example we consider the deformation of the unit square [0, 1] × [0, 1] as in Sec.1.1.3( see Fig. 8).The set of varying parameters ρ is the same as in Sec.1.1.3.The reference domain and the associated mesh T ref , the generic domain and associated mesh T h are also the same as in Sec.1.1.3.The functional F ρ used to describe the boundary Γ ρ is defined by: In Fig. 12 we represented the reference mesh (left) T ref and the deformed mesh T map with α = 2 and β = 1 (right).As previously, we observe that the deformed mesh fits well with the objective in green.
We consider the same Laplace problem as in the example 1.1.3.In Fig. 13 we represented the solution of the Laplace problem φ map (left) and the relative error measured in the L ∞ -norm between φ h and I h φ map , where I h is the interpolation operator from X map into X h .As expected, this error is on the order of h 2 .ρ and Γ 3 ρ is defined similarly as in Sec.1.2.1.However, now on the boundary Γ 2 ρ the functional is nonlinear (see Fig. 14) and defined by: In addition to the coefficients α and β, the amplitude will also belong to the set of varying parameters ρ.

Γ3
ref Because of the nonlinearity of F ρ , we use a Picard fixed-point algorithm to solve problem (18) that can be rewritten under the form or equivalently, defining the solution operator Starting from an initial guess u u u 0 , we solved iteratively the following problem for n = 1, In Fig. 15 we show the reference mesh T ref (left), and the deformed mesh T map for α = 1, β = 0.7, and = 0.1 (right).We observe that the deformed mesh fits well with the objective in green.A better quantification of the quality of the meshes, using quantities associated to the mesh as in Sec.1.1.1,is given in Table 3.These results highlight the regularity of the mesh obtained by our transformation, for a non-uniform shear square (with nonlinear F ρ ).In this example, we consider a Laplace problem (1) where g f is set to g f (x x x) = x 2 (1 − x 2 ) on Γ f .In Fig. 16 we show the solution of the Laplace problem φ map (left) and the relative error measured in the L ∞ -norm between φ h and I h φ map .As previously, this error is on the order of h 2 (see Table 3), as in the previous examples.r ho is also defined similarly as in Sec.1.2.2.However, now on the boundary Γ 2 ρ the functional F ρ is defined by: , where the coefficient α represents the set of varying parameters ρ (see Fig. 17).

Γ2
ref To treat the nonlinearity of F ρ , we choose to use a steepest-descent method with step size ζ > 0 to solve problem (18).Starting from an initial guess u u u 0 , we iteratively compute for n = 0, 1, • • • , n max : where d d d n is the solution to the problem: Find Such approach may lead to inverted triangles during the process, thus we propose the following alternatives for the term (d d d n , v v v) elas which improve the method: where φ(t) = 1 − e −t , and κ is the aspect ratio of triangles, defined by κ = rcirc 2 * rin , where r circ is the radius of the circumscribed circle, and r in is the radius of the incircle of the triangle.
This idea of changing the inner product whereby a gradient is identified for J(w w w) is quite classical in the study of gradient flows; it amounts to an efficient preconditioning of the minimization problem (16); see e.g.[4] about this point.
In Fig. 18 we represented the reference mesh T ref (left), and the deformed meshes for α = 0.3, using (., .)elas,F (middle) and (., .)elas,κ (right).We observe that the mapped mesh fits well with the objective for both cases and that the triangulation obtained with (., .)elas,κ is more regular than the one obtained with (., .)elas,F .

Deformation of a crescent moon
In this example, we consider a generic component in the form of a square [−1, 1] × [−1, 1] where a crescent moon has been removed.The crescent moon is the intersection of the exterior of a disc of radius r 1 centered in (c, 0.5) with a disc of radius r 2 centered in (0.5, 0.5) (see Fig. 19), where the coefficients c, r 1 and r 2 represent the set of varying parameters ρ.In the generic domain Ω ρ , the parameter dependent boundary Γ ρ is made of the union of Γ (1, 1) The functional that describes the boundary Γ ρ is nonlinear and as follows: To treat the nonlinearity of F ρ in (18) we choose to use one of previous steepest descent algorithm with (., .

Recap. on the Reduced Basis Method
In this section we present the mapping used in a reduced basis method context.The implementation is then decomposed in the following independent steps.

First step: construction of the manifold of deformations
Following one of the methods presented above, we are able to solve an elasticity problem, the solution of which is a displacement that enables to go from Ω ref onto the domain of interest Ω ρ and that maps the points of each part Γ m ref of the boundary of Ω ref into (actually onto) the associated part Γ m ρ of Ω ρ .We compute such elasticity solutions for a large number N of values of ρ hopefully representing well the set of all problems we shall be faced to (in our case N was set to 100).The displacements are denoted as U U U (ρ), these are defined over Ω ref , the mapping from It is expected, verified in our applications -and it would be good to prove it -that the set of all {u u u m (ρ)}, when ρ varies is a manifold with a small Kolmogorov n-width, which is the requirement for next building a sensible reduced basis approach (see e.g.[9,17]).

Second step: extraction of a reduced basis for fast approximation of the deformations for general parameters
We extract with a POD or a greedy procedure (see e.g.[9,17]) from this manifold {u u u m (ρ)} when ρ varies, a (small) set of parameters ρ 1 , ρ 2 , . . ., ρ n , . . .such that, for any given ε > 0, there exists n = n(ε) such that, for any ρ, there exists components α 1 (ρ), α 2 (ρ), . . ., α n (ρ) such that u u u m (ρ) − where C is some stability constant.Actually, the coefficients α i (ρ) can be found by many ways from the knowledge of the N solutions u u u m (ρ) that were computed, but it is also possible to get them for parameters that do not belong to the set of parameters that have been chosen in subsection 2.1.If these chosen values indeed represent well the set of all problems we shall be faced to, then we can propose to define, for any ρ the {α i (ρ)} i=1,...,n by least square applied to the implicit caracterization F ρ (.) = 0 defining the boundary of Ω ρ .This means that {α i (ρ)} i=1,...,n = arg min {βi}i=1,...,n F ρ n i=1 β i u u u m (ρ i ) L 2 (∂Ω ref ) (20) Having found these values {α i (ρ)} i=1,...,n , an approximation of the elasticity problem that maps Ω ref over Ω ρ is thus also given by n i=1 α i (ρ) U U U (ρ i ).The way we solve (19) and ( 20) is as follows : for (19), an EIM Greedy approach is used to identify the u u u and hence the U U U .For (20), we used a standard matlab nonlinear least-squares solver ('lsqnonlin').Note that (20) produces the optimal coefficients, and is thus not a POD approach that targets to find the optimal u m (ρ i ).The Table 4 summarize the computational times and error of the reduced basis method.Here the errors are absolute.However, the solutions are O(1) so relative and absolute are quite similar.

Case of a bell
We consider a 2D test case of a bell.The geometrical parameter is α (see Fig. 17).In Fig. 22 we represented the solution of the Laplace problem (1) by coupling the mapping method to the reduced basis method.In Table 5, the computational times and error of the reduced basis method are shown.Here the errors are absolute.However, the solutions are O(1) so relative and absolute are quite similar.

Figure 2 :
Figure 2: (1+ν) ) are the Lamé coefficients, with E the Young's modulus and ν the Poisson's ratio, and e e e is the linearized strain tensor given by e e e 11 (u u u) = ∂ x1 u 1 , e e e 22 (u u u) = ∂ x2 u 2 , e e e 12 (u u u) = e e e 21 (u u

Figure 5 :
Figure 5: Left: solution φ map ; Right: relative error between φ h and I h φ map .

Figure 7 :
Figure 7: Left: solution φ map ; Right: relative error between φ h and I h φ map .

Figure 9 :
Figure 9: Left: reference mesh T ref ; Right: deformed mesh T map .

Figure 11 :
Figure 11: Left: solution φ map ; Right: relative error between φ h and I h φ map .

Figure 12 :
Figure 12: Left: reference mesh T ref ; Right: deformed mesh T map .

Figure 13 :
Figure 13: Left: solution φ map ; Right: relative error between φ h and I h φ map .

Figure 15 :
Figure 15: Left: reference mesh T ref ; Right: deformed mesh T map .

Figure 16 :
Figure 16: Left: solution φ map ; Right: relative error between φ h and I h φ map .
with (w w w, v v v) elas := Ω ref 2µe e e(w w w) : e e e(v v v) + λdiv (w w w)div (v v v)dx x x.
) 2µe e e(d d d n ) : e e e(v v v) + λdiv (d d d n )div (v v v) dx x x (d d d n , v v v) elas,κ := Ω ref 1 φ( 1 κ γ ) 2µe e e(d d d n ) : e e e(v v v) + λdiv (d d d n )div (v v v) dx x x, with γ = 7,

α
i (ρ) u u u m (ρ i ) L ∞ (∂Ω ref ) ≤ ε (19)By linearity of the elasticity problem, the solution U U U (ρ) to the elasticity problem over Ω ref , with Dirichlet boundary conditions u u u m (ρ) is thus close to n i=1 α i (ρ) U U U (ρ i ) with an error over Ω ref that is bounded by Cε

Figure 21 :
Figure 21: Solution of Laplace problem

1
Computation of the displacement field for the mapping and resolution of the Laplace equation on the deformed domain Let ρ be a set of geometrical factors used to parameterize the geometry of the unique component C 1,ρ ≡ Ω ρ and Ω ref the reference domain.The aim of this paper is to propose automatized techniques to build a mapping T T T from the domain Ω ref to the domain Ω ρ .We consider here only the two-dimensional case (d = 2): we denote by (x 1 , x2 ) the coordinates of the point x x x in the reference domain Ω ref , (x 1 , x 2

Table 1 :
Example 1: classical quantities associated to the meshes T ref , T map and T h

Table 2 :
1 and Table2are almost the same, and thus using penalized conditions appears to be a good alternative to using strong boundary conditions.Example 2: classical quantities associated to the mesh T map on Γ m ref can be considered as a Dirichlet boundary condition for an elasticity problem that maps Ω ref onto the domain of interest Ω ρ .In opposition to what generally happens for those Dirichlet boundary conditions, they are not imposed a priori but are obtained a posteriori, after the problem has been solved.