Luca Cirrottola

Mario Ricchiuto

Algiane Froehly

Barbara Re

Alberto Guardone

Giuseppe Quaranta

Adaptive deformation of 3D unstructured meshes with curved body fitted boundaries with application to unsteady compressible flows

Keywords: Constant-connectivity mesh adaptation, unstructured meshes, unsteady compressible flows, conservative formulations

We present an adaptive moving mesh method for unstructured meshes which is a threedimensional extension of the previous works of Ceniceros et al.

Ceniceros2001

[10], Tang et al.

Tang2003

[40] and Chen et al.

. The iterative solution of a variable di↵usion Laplacian model on the reference domain is used to adapt the mesh to moving sharp solution fronts while imposing slip conditions for the displacements on curved boundary surfaces. To this aim, we present an approach to project the nodes on a given curved geometry, as well as an a-posteriori limiter for the nodal displacements developed to guarantee the validity of the adapted mesh also over non-convex curved boundaries with singularities.

We validate the method on analytical test cases, and we show its application to two and three-dimensional unsteady compressible flows by coupling it to a second order conservative Arbitrary Lagrangian-Eulerian flow solver.

Introduction sec:Intro

Mesh adaptation is a powerful tool to improve the representation of complex fields for a given computational expense. In computational fluid dynamics in particular, adaptation has become nowadays a customary tool Park2016 [START_REF] Park | Unstructured grid adaptation: Status, potential impacts, and recommended investments towards cfd 2030[END_REF]. Adapting the mesh also has a relative computational overhead, which motivates the quest for e cient and robust methods.

Techniques improving the discrete representation of the fields of interest by inserting and removing mesh entities (so called h-adaptation) have proven to be quite mature Park2016 [START_REF] Park | Unstructured grid adaptation: Status, potential impacts, and recommended investments towards cfd 2030[END_REF].

However, solution transfer between meshes with di↵erent topologies may be non-trivial and may have a non-negligible computational cost, especially if conservation constraints need to be satisfied FarrellMaddison2011,Alauzet2016,hermes2018highorder,KucharikShashkov2008,Re2017 [START_REF] Farrell | Conservative interpolation between volume meshes by local galerkin projection[END_REF][START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF][START_REF] Hermes | High-order solution transfer between curved triangular meshes[END_REF][START_REF] Kucharik | Extension of e cient, swept-integration-based conservative remapping method for meshes with changing connectivity[END_REF][START_REF] Re | An interpolation-free ale scheme for unsteady inviscid flows computations with large boundary displacements over three-dimensional adaptive grids[END_REF].

By constrast, mesh nodes relocation with constant element connectivity (so called r-adaptation), o↵ers the possibility of a minimally intrusive coupling with existing com-putational mechanics solvers, as no modification of the data structures is required. As h-adaptation methods, they also provide considerable improvement in the quality of the solutions obtained, especially in unsteady simulations of traveling waves, like shock waves and water waves, where uniform refinement would be way too costly. Moreover, with r-adaptation methods devising conservative projections is much simpler. In fact, the preservation of the one-to-one mapping from the old to the new mesh entities allows the easy construction of a conservative remapping Kucharik2003efficient [START_REF] Kucharik | An e cient linearity-and-boundpreserving remapping method[END_REF], or the use of Arbitrary-Lagrangian-Eulerian (ALE) formulations compliant with the Geometric Conservation Law (GCL)

Trulio1961,ThomasLombard1979 [START_REF] John | Numerical solution of the one-dimensional lagrangian hydrodynamic equations[END_REF][START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF].

Unfortunately, the preservation of the initial mesh topology undeniably imposes severe restrictions on nodal displacements in order to avoid mesh folding with tangled (i.e. inverted) elements, especially when the boundary exhibits singular points. Moreover, the accuracy attainable for complex solutions is limited by the initial density of mesh nodes, and less finely-tunable than in metric-based h-adaptation AlauzetFrey2003 [START_REF] Alauzet | Estimateur d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage. Partie I : aspects théoriques[END_REF].

Anyway, the advantages brought by the e↵ortless coupling with external flow solvers and the conservative solution remapping can counterbalance the mesh quality limitation as long as the r-adaptation technique is computationally e cient. Simply put, the error reduction brought by adapting the mesh must o↵set the computational overhead. A measure of this e ciency can be evaluated by comparing with the cost of a simulation run on a uniformly refined mesh, providing the same resolution of the flow field. Hybrid adaptive approaches combining well timed re-meshing and adaptive deformation at every time step, which are perhaps the ones computationally most appealing, still require the r-adaptation step to perform well.

Extensive reviews of r-adaptation can be found in Knupp1993,Budd2009 [START_REF] Knupp | Fundamentals of grid generation, The Fundamentals of Grid Generation[END_REF][START_REF] Budd | Adaptivity with moving grids[END_REF]. We focus here on methods based on the numerical solution of an elliptic partial di↵erential equation for the position of the mesh nodes, often referred to as the mesh PDE. This equation is typically formulated to find a mapping ⇠ : ⌦ x ! ⌦ ⇠ from the physical domain to a reference (computational) one. This mapping needs to be injective and surjective in order to guarantee that the produced mesh neither folds nor breaks the domain. Historically the Winslow or homogeneous Thompson-Thames-Mastin generator TTM1977,Thompson1977boundary [START_REF] Joe F Thompson | Tomcat -a code for numerical generation of boundary-fitted curvilinear coordinate systems on fields containing any number of arbitrary two-dimensional bodies[END_REF][START_REF] Joe F Thompson | Boundary-fitted curvilinear coordinate systems for solution of partial di↵erential equations on fields containing any number of arbitrary two-dimensional bodies[END_REF] x ⇠ = 0 has been the basis for structured boundary-fitted grid generation (see for example the review in Thompson1985 [START_REF] Thompson | Numerical grid generation: foundations and applications[END_REF])

and it has been extended to also adapt the mesh in the domain either by adding source terms to the equation, or through a variable-di↵usion approach Winslow1981 [START_REF] Winslow | Adaptive-mesh zoning by the equipotential method[END_REF]. A more general formulation of the last method has been given in Dvinsky1991 [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF] by means of harmonic maps and extended in Brackbill1993 [START_REF] Brackbill | An adaptive grid with directional control[END_REF].

These equations describe the mapping ⇠ = ⇠(x) and need to be inverted for the physical coordinates x = x(⇠), leading to a nonlinear system of PDEs which is iteratively solved. In order to ease the cost of the iterative solution of the inverted equations, an alternative mesh generator was proposed in Ceniceros2001 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF] based on a variable-di↵usion Laplace equations directly formulated in the physical domain for the mapping x : ⌦ ⇠ ! ⌦ x . This generator is not based on a theoretical derivation, but on the observation that the variable-di↵usion Laplacian in the reference domain is su cient to adapt the mesh in the desired regions while the equations, which are still nonlinear, can easily be solved through a relaxation procedure. The e ciency of the method was shown in application to two-dimensional Boussinesq convection on structured grids, and the method was later applied in Tang2003 [START_REF] Tang | Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws[END_REF] to hyperbolic conservation laws and extended in Chen2008 [START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF] to multicomponent flows on two-dimensional unstructured grids. More recently, the same method was applied to the two-dimensional shallow water equations both in Cartesian and spherical coordinates arpaia:hal-01102124,arpaia:hal-01372496,Arpaia2020 [START_REF] Arpaia | Mesh adaptation by continuous deformation. basics: accuracy, e ciency, well balancedness[END_REF][START_REF] Arpaia | r-adaptation for shallow water flows: conservation, well balancedness, e ciency[END_REF][START_REF]Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF].

Robustness to mesh folding in r-adaptation is a delicate matter, similarly critical in the context of mesh deformation related to moving boundaries, curved mesh generation and smoothing (see for example JohnsonTezduyar1994,Dwight2009,Toulorge2013,Fortunato2016,Turner2018 [START_REF] Johnson | Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces[END_REF][START_REF] Dwight | Robust mesh deformation using the linear elasticity equations[END_REF][START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF][START_REF] Fortunato | High-order unstructured curved mesh generation using the winslow equations[END_REF][START_REF] Turner | Curvilinear mesh generation using a variational framework[END_REF]). Obtaining non-singular meshes requires two main conditions to be met. The first is that the continuous map, appropriately modified to account for all boundary conditions, should verify the appropriate conditions as e.g. the non-negativity of the determinant of the deformation Jacobian.

Until quite recently, su cient conditions were known only in the framework of harmonic maps Dvinsky1991,Liao1992 [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF][START_REF] Liao | Variational approach to grid generation[END_REF]. Recent work by Huang2018 [START_REF] Huang | On the mesh nonsingularity of the moving mesh pde method[END_REF], has allowed to prove similar properties for other types of mesh PDEs, as e.g. some of those proposed by Huang Huang2001 [START_REF] Huang | Variational mesh adaptation: Isotropy and equidistribution[END_REF] or Huang and Russel Huang2011 [START_REF] Huang | Adaptive moving mesh methods[END_REF], by resorting to energy arguments borrowed from the theory of gradient flows.

The second important aspect is that the discretization used to approximate the mesh PDE should have the appropriate "property preserving" character, so that the fully discrete moving mesh method is also guaranteed to provide non-singular meshes. This is in itself a subject of investigation. It is in general well known that discrete moving mesh methods can lead to mesh tangling even with properly chosen mesh PDEs Dvinsky1991 [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF], and the impact of the truncation error is stressed for example in Knupp1993 [START_REF] Knupp | Fundamentals of grid generation, The Fundamentals of Grid Generation[END_REF]. In the setting of gradient flow maps, geometrical discretizations have been shown in Huang2018 [START_REF] Huang | On the mesh nonsingularity of the moving mesh pde method[END_REF] to answer the discrete positive Jacobian requirement. However, even in the last reference, the issue of accounting for complex curved boundaries is overlooked, even though mesh movement along a given surface does not appear to be necessarily a natural boundary condition of the PDEs considered.

In this work we proceed di↵erently. We want to be able to handle domains with boundaries as general as possible in three space dimensions. and propose a relaxation technique embedding a geometrical limiter allowing to achieve this objective. We focus on the simple reference-domain variable-di↵usion Laplacian approach originally proposed in Ceniceros2001 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF], however the ideas proposed in this paper can be extended to other mesh PDEs. To the best of the authors' knowledge, very few applications of r-adaptation to three-dimensional meshes are available to date, see for example Li2002 [START_REF] Li | A moving mesh finite element algorithm for singular problems in two and three space dimensions[END_REF] as well as some simple applications in Huang2018 [START_REF] Huang | On the mesh nonsingularity of the moving mesh pde method[END_REF]. The original approach in Ceniceros2001 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF] does not o↵er theoretical guarantees against folding, although it has been successfully used in many applications with non-convex boundaries. R-adaptation in three dimensions exhibits even stronger limitations than in two dimensions, as su ciency results for unfolded continuum maps are typically based on requirement of smoothness and convexity of the boundary. These are easily violated especially in application to external flows, where the boundary is not convex and several singularities (like corners and ridges) are present. We have experienced that tangling is a major concern in the three-dimensional extension of these techniques.

Smoothing or untangling methods for unstructured meshes already developed in the literature require nontrivial procedures Hansen2004,Toulorge2013 [START_REF] Hansen | A finite element method for unstructured grid smoothing[END_REF][START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF].

Our contribution is thus related to a mesh displacement method allowing to guarantee that nodes move and always remain on a given parametrizaion of curved domain boundaries, and for an a-posteriori limiter for the nodal displacement which, when embedded in the mesh relaxation iterations, allows to prevent the occurrence of tangled elements, thus enforcing the validity of the discrete mapping while avoiding smoothing procedures. The resulting moving mesh library Fmg has been developed on top of the open source software Mmg Dapogny2014,MmgSW [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF]Mmg software for simplicial remeshing[END_REF] to exploit, among other things, its built-in cubic Bézier patch representation of complex manifolds.

The paper is organized as follows. We recall the continuous mesh partial di↵erential equations in section §2, while a thorough discussion of their numerical solution is given in section §3, including a-posteriori limiting and projections to obtain a valid mesh satisfying all the boundary conditions, and the application to unsteady simulations.

We discuss the validation of the method proposed considering the adaptation w.r.t. analytical functions in two and three space dimensions in section §4, and application to two and three dimensional unsteady compressible flows are discussed in section §5 Finally, conclusions are presented in section sec:Conclusions 6.

Variable-di↵usion Laplacian r-adaptation in the reference domain sec:NumModel

We focus here on Laplacian-based r-adaptation, which is the mesh PDE currently implemented in the Fmg library. However, the ideas proposed in this paper can be immediately extended to other mesh PDEs. We recall here the continuous mesh problem, and in particular we discuss the boundary conditions, as well as the definition of the monitor functions used for adaptation.

Following

Chen2008 [START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF], we look for a mapping x : ⌦ ⇠ ! ⌦ x from the reference domain ⌦ ⇠ (the original mesh) to the computational domain ⌦ x (the adapted mesh). Within the reference domain, the computational coordinates satisfy the variable-di↵usion Laplace equation

r ⇠ • (!(x)r ⇠ x) = 0 in ⌦ ⇠ (1) eq:rLap
The above problem is in general a system of coupled non-linear PDEs, which needs to be complemented by appropriate boundary conditions. Nonlinearity is introduced by the monitor function !(x) which depends on an external field evaluated in the computational domain. In this work, we have used a classical scalar definition for this quantity, which allows to decouple the equations for the three spatial coordinates.

In particular, given a quantity of interest f (x), the scalar monitor function used in the examples discussed later is evaluated as

!(x) = q 1 + ↵||r ⇠ f (x)|| 2 ↵ + ||H ⇠ (f )(x)|| 2 + ⌧ ||f || 2 ⌧
(2) eq:monitorFunction where r ⇠ and H ⇠ denote the gradient and Hessian computed on the reference domain

⌦ ⇠ . Norm ||f || is defined as ||f || = min ✓ 1, ||f || max(||f ||) ◆ . (3)
This normalization, already used in Chen2008 [START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF], allows to introduce some saturation near the norm maximum according to the value of . The idea behind this normalization is to spread a little bit the peak values of the function f around the peak locations, in order to filter out small inhomogeinities in the numerical approximation of the sharp fronts of f . The above definition gives the user some control on the behaviour of the spatial mapping via the parameter pairs (↵, ↵ ), ( , ), (⌧, ⌧ ). As in the original works Ceniceros2001,Chen2008 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF][START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF],

these parameter pairs are not related to an error estimate, thus they are determined empirically and they are intrinsically dependent on the application. However, in our experience a short test on a few time steps is su cient to assess the behaviour of these parameters on the whole simulation time range.

Boundary conditions sec:NumModelBCs

Despite the decoupling of equation ( eq:rLap 1) into separate scalar equations, even for a scalar monitor function a strong coupling of the coordinate equations may arise through the boundary conditions, especially in domains including general shapes. In particular, we will split the boundary on two parts as

@⌦ ⇠ = D ⇠ [ S ⇠ . A full set of Dirichlet conditions are imposed on D ⇠ x = ⇠ on D ⇠ ( 4
)
as this is the portion of the boundary that is not allowed to move. Since fixed boundary nodes are not convenient whenever adaptation is performed on flow waves approaching the boundary, we want to limit as much as possible usage of Dirichlet conditions, preferring slip boundary conditions wherever we can formulate them while preserving sharp geometrical features of the boundary. Along the slip boundary S ⇠ the coordinate positions are constrained to move along a given parameterized domain. For manifold surfaces, we assume to have a known parameterization, for example in the form S (x) = 0.

This provides one position constraint relating the d spatial coordinates. Thus, d 1 additional conditions are required, which are here taken as the null normal stress conditions parallel to the local tangent space spanned by {⌧ S j } j=1,d 1 . This gives the boundary conditions:

S (x) = 0 nS • (!(x)rx) • ⌧ S j = 0 8j = 1, d 1 on S ⇠
(5) eq:manifold-bc

In three dimensions, the slip boundary S ⇠ can be further generalized as the union of multiple manifolds with one-dimensional boundaries joining them. This leads to intersection curves where two sets of equations eq:manifold-bc 5 should formally be satisfied at the same time, constraining the displacement to happen along the curve tangent direction.

Points where multiple intersection curves meet are corner points and no displacement is possible, thus Dirichlet conditions are imposed on (and only on) these points. The approach used here for the numerical enforcement of boundary conditions is discussed in the following section.

3. Discrete equations, a-posteriori limiting, slip on curved boundaries sec:NumModel1

We discuss here the implementation choices made in the Fmg library, namely the discretization of the mesh PDEs, as well as their iterative solution. Both the a-posteriori limiting of the displacement and the implementation of the slip boundary conditions are strongly tied to the relaxation iterations, and for this reason all the steps are discussed in this section. More specifically, to relieve the complexity of satisfying the boundary conditions ( eq:manifold-bc 5), we formulate the discrete approximation by means of an iterative multiple-corrections procedure embedding the following three elements:

1. A finite element approximation of the variational form of ( eq:rLap 1) with natural (Neumann) boundary conditions;

2. An a-posteriori limiter for the nodal displacement enforcing local mesh valididy;

3. A boundary correction in the local normal direction to enforce the first of ( eq:manifold-bc

5) by

projecting on the manifold parametrization at hand.

The intertwining of the a-posteriori limiter of the displacement, of the update of the local boundary normals, and of the projection on the parameterized manifold is essential for the proposed approach to provide valid adapted meshes both in the volume and on the boundaries.

Concerning the representation of these, the Fmg library we developed makes use of the point-normal curved triangles parameterization proposed in Vlachos2001 [START_REF] Vlachos | Curved pn triangles[END_REF], which is based on cubic Bézier patches with quadratically varying normals. For this we leverage the implementation provided by the open source software Mmg Dapogny2014,MmgSW [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF]Mmg software for simplicial remeshing[END_REF]. However note that the form of the manifold parametrization is not a necessary ingredient of our method.

Other high order approximations can be used.

Finite element approximation: Dirichlet and natural boundary conditions sec:weak

The discrete equations are built starting from a linear finite element approximation of the problem embedding natural (Neumann) boundary conditions, which corresponds to the simple variational form

Z ⌦ ⇠ !(x)r ⇠ v • r ⇠ x d⌦ ⇠ = 0, 8v 2 H 1 (⌦ ⇠ ) . ( 6
) eq:var-rlap
Note that this statement satisfies the null normal stress conditions in ( eq:manifold-bc 5) but neither the remaining one in the system (the belonging to the surface), nor any Dirichlet conditions eventually assigned. Dirichlet conditions are strongly imposed on the solution space, while the enforcement of the belonging to a slip surface will be detailed in the following sec:slip 3.4) Note also that ! being a scalar quantity, the above equations provide uncoupled nonlinear variational statements for each component of x.

The projection of ( eq:var-rlap 6) on the linear finite element space leads to the nonlinear algebraic system [START_REF] Brackbill | An adaptive grid with directional control[END_REF] eq:fem1

K(x)x ⌫ = 0, ⌫ = 1, . . . , d (
having introduced the array of unknown node positions x ⌫ = [x ⌫ i ] for each space component ⌫, with d the number of space dimensions, and where the sti↵ness matrix has the standard entries

K ij (x) = Z ⌦ ⇠ !(x)r ⇠ i • r ⇠ j d⌦ ⇠ (8) eq:stiffLap
with { i } i 1 the linear base functions spanning the solution space. Please note that ( eq:fem1 7)

is a set of decoupled systems, one for each spatial direction, as shown by the fact that K ij are scalar entries. Note also due to the consistency of the finite element space with Dirichlet conditions, Dirichlet nodes are not included in the above sum.

In practice, by defining the displacement = x ⇠, the system is not written as in

( eq:fem1 7
), but as

K(x) ⌫ = K(x)⇠ ⌫ (9) eq:fem2
which is better suited for the iterative corrections described in the following sections.

Scalar correction iterations sec:Jacobi

We introduce an iterative procedure which, while avoiding mesh tangling in 2D and 3D, and accounting for the directional coupling inherent to ( eq:manifold-bc 5), retains the scalar structure of the decoupled variational form. Note however that the corrections proposed can be easily adapted to other iterative solution methods (as well as mesh PDEs).

The basic iteration used in our method starts from a standard diagonal Jacobi relaxation to handle the nonlinearity of ( eq:fem2

9

)

K [k] ii [k+1] i = X j2B i j6 =i K [k] ij [k] j X j2B i K [k] ij ⇠ j ( 10
)
where B i denotes the ball 1 of node i and vector

[k] i = [ ⌫ i ] [k]
⌫=1,...,d is now the vector made of the space components of the displacement of node i at iteration k (the same notation will be used for vectors x i and ⇠ i ). Again, we stress that the above iteration is in fact a set of d relations for the components of the displacement. The matrix entries K

[k] ij depend on monitor function ! at iteration k (cf. equation (

eq:stiffLap 8
)), which in turn depends on the scalar field f evaluated at the actual positions x In our implementation we added and removed the term

K [k] ii [k]
i , to obtain the following iterations

[k+1] i = [k] i 1 K [k] ii X j2B i K [k] ij x [k] j ( 11
)
1 In the common lexicon of the mesh generation community, the ball Bi of a mesh node i is defined here as the set of elements sharing vertex i FreyGeorge2007 [START_REF] Pascal | Mesh generation: Application to finite elements[END_REF]. To ease the notation, with a slight abuse of terminology we say that a node j belongs to the ball of node i if it is a vertex of an element belonging to the ball, and we denote this statement as j 2 Bi.

which are initialized with

[0] i = 0.
The last step is the computation of the new nodal positions as follows

x [k+1] i = x [k] i + f x [k+1] i ⇣ [k+1] i , {x [k+1] i } j<i , {x [k] i } j i ⌘ ( 12
)
where

g x i [k+1]
are limited increments obtained by a-posteriori correcting

[k+1] i
to account for both mesh validity and boundary conditions (both Dirichlet and slip wall).

In both cases, these corrections are local, albeit not only dependent on node i, and non-linear w.r.t. x. The nonlinearity is readily handled in the iterations by using the last nodal positions available.

A-posteriori corrections for mesh validity enforcement sec:relax

The Laplacian model in the reference domain does not guarantee that the Jacobian of the mapping is strictly positive everywhere, thus leading to the occurrence of tangled (invalid) mesh elements. In two space dimensions, our experience has shown that in most cases carefully tuning the monitor function !(x) allows to solve this issue. This is not the case in three dimensions, where tangling occurs much more often.

To cope with this, we have devised an a-posteriori limiter to the nodal displacements which is activated whenever the displacement of a node causes the occurrence of an element whose volume is below a given threshold. This condition is of course implicit, in the sense that it couples the positions of all the mesh nodes. However, it can be easily embedded in an iterative setting. In particular, in our implementation we relax and update each nodal position x

d 0 i = [k+1] i + ⇠ i x [k] i d s+1 i = ( µ i d s i if min K2B [k,k+1] i |⌦ K | < ✏ d s i otherwise 8s 2 [0, . . . , s max 1] f x [k+1] i = d smax i ( 13
) eq:limiter
The limiter is thus the result of local sub-iterations, which are stopped when a volume greater than ✏ is guaranteed for every element in the ball. This check allows to enforce the validity of every intermediate mesh configuration, e↵ectively preventing the occurrence of invalid elements at a reasonable computational cost with respect to smoothing or untangling procedures for unstructured meshes

Hansen2004,Toulorge2013 [START_REF] Hansen | A finite element method for unstructured grid smoothing[END_REF][START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF].

It must be remarked that for interior nodes in general one iteration of the above procedure is enough, while more iterations are required when applying the limiter within the projection step enforcing the boundary conditions (cf. next section). For simplicity here the same value µ i = 0.5 has been adopted for all the nodes. In the present work, this value has proven to be e↵ective in locally preventing tangling while allowing the position of blocked nodes to be naturally relaxed by the successive application of r-adaptation at the next time steps, as the monitor function moves with the flow.

x The decoupling of the spatial coordinates obtained by initially accounting for natural boundary conditions only is particularly convenient in terms of computational cost and simplicity of implementation. It allows to store and assembly only a single smaller sti↵ness matrix to be used for every space coordinate, instead of a matrix of 3 ⇥ 3 blocks. However, the resulting nodal displacements need to be corrected to account for conditions on Dirichlet and slip boundaries. This is achieved by the projection step discussed in this section, which is easily embedded in the scalar iterations. The description is given for slip wall boundaries, of which Dirichlet nodes are a particular case.

[k] i x [k] i + d 0 i fig:relaxStart (a) Invalid displacement. x [k] i x [k] i + d 1 i fig:relaxCorrect (b) Relaxed displacement. x [k+1] i fig:relaxOk (c) Updated configuration.
In the Fmg library boundary geometries are handled by means of curved point-normal triangles Vlachos2001

[51], i.e. piecewise cubic Bézier patches for the boundary position and quadratic for the boundary normal vector, relying on the implementation provided in Mmg Dapogny2014,MmgSW [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF]Mmg software for simplicial remeshing[END_REF].

In this setting an implicit surface representation of the slip boundary reading in the continuous case [START_REF] Drikakis | Computation of non-stationary shockwave/cylinder interaction using adaptive-grid methods[END_REF] eq:surface is approximated by the explicit piecewise parametric representation [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF] eq:model_geo which is defined for each triangle ⌧ in the triangulation of the slip boundary, from the positions and unit normals {x i , ni } i2⌧ of the nodes of the triangle. Similarly, the Bézier patches also allow to evaluate surface normals as

S (x) = 0 on S ⇠ (
⌧ [x j , nj ] : ⌃ ! S ⇠ , ⌃ = [0, 1] ⇥ [0, 1], S ⇠ ⇢ R 3 x = ⌧ [x j , nj ](w), w 2 ⌃ (
⌘ ⌧ [x j , nj ] : ⌃ ! S ⇠ , ⌃ = [0, 1] ⇥ [0, 1], S ⇠ ⇢ R 3 n = ⌘ ⌧ [x j , nj ](w), w 2 ⌃ (16) eq:model_normal
For some applications, as for example external aerodynamics, handling curved geometries is a necessity. As a consequence, the geometric approximation becomes an integral part of the numerical method. In particular, in three dimensions even the simplest combination of boundary surfaces easily leads to intersection curves. Since sharp edges (ridges) in the initial geometry need to be preserved as well as corners (intersections of multiple ridges), nodes cannot cross a ridge, but they are only allowed to move tangentially to it, and displacement of a corner node cannot happen. As outlined in sec- , which allows to build a local updated geometrical model. As before, this model is evaluated using the new updates for nodes already processed, and value from the previous iteration for the remaining ones. This provides the incrementally updated geometry model

[k,k+1] ⌧ = ⌧ [{x [k+1] j , n[k+1] j } j<i , {x [k] j , n[k] j } j i ] (cf. ( eq:model_geo 15
)). In particular, as shown on figures The iterations providing the final displacement, and hence position, are similar to ( eq:limiter 13):

d 0 i = [k,k+1] ⌧ ⇣ w(x [k+1] i ) ⌧ ⌘ x [k] i d s+1 i = 8 < : [k,k+1] ⌧ ⇣ w(x [k] i + µ i d s i ) ⌘ x [k] i if min K2B [k,k+1] i |⌦ K | < ✏ d s i otherwise , 8s 2 [0, . . . , s max 1] f x [k+1] i = d smax i (17) eq:relaxSlip
We stress again that since the piecewise patches depend on both node positions and unit normals, the position update is always accompanied by the re-evaluation of the unit normal vectors through the analogously defined model ⌘

[k,k+1] ⌧ (cf. ( eq:model_normal 16
)). This is omitted from ( eq:relaxSlip [START_REF] Ashley | An evaluation of several di↵erencing methods for inviscid fluid flow problems[END_REF]) to keep a lighter notation.

Ridges. The displacement check and projection on boundary ridges is handled exactly in the same way as for manifold surfaces. The main di↵erence is that now the parametric space is replaced with a curve parametrisation which is one-dimensional ⌃ ⇢ R. Thus all operations previously performed on the tangent plane are performed by projection on the tangent line, and normal vectors of both the manifold surfaces joining at the ridge are stored and updated in the geometrical model.

Corners. These are the only allowed Dirichlet nodes, thus corners verify exactly the boundary condition, and are not included in the discrete variational form sec:weak 3.1. In this specific case, displacement is not allowed as they are already on the exact geometry, and the condition imposed is

x [k+1] i = ⇠ i (18)
3.5. Unsteady mesh adaptation through restarted iterations eq:r:adaptDyn

Following

Tang2005,Chen2008 [START_REF] Tang | Moving mesh methods for computational fluid dynamics[END_REF][START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF], dynamic mesh adaptation during the time evolution of a fluid flow simulation is performed by repeating the steady adaptation procedure described in the previous section at each time step, without the explicit formulation of a di↵erential equation in time for mesh motion. This simplifies the coupling with existing flow solvers.

In this case of fixed boundary domains, the reference mesh ⇠ is constant in time, while the computational mesh x(t (n+1) ) is the r-adaptation of the (fixed) reference mesh.

Thus, the displacement at each time step n + 1 is initialized with the value achieved at the last Jacobi iteration K achieved in the previous time step n

[0](n+1) i = [K](n) i (19)
so that successive Jacobi iterations during time evolution are e↵ectively accumulated on the nodes positions.

Validation via adaptation on analytical functions sec:AnalyticalRes

We consider here a series of analytical tests allowing to measure the e↵ectiveness of the method. As shown in section sec:NumModel 2, we recall here that the mesh adaptation model can be governed by the number of iterations n it plus the three parameter pairs (↵, ↵ ), ( , ), (⌧, ⌧ ), representing the intensity and the normalization constant of the solution gradient, the solution Hessian, and the solution itself in the definition of the monitor function. In this work, we have not seen specific benefits in mixing all three parameter pairs, so we will explicitly report only the values for the used pairs, while values not shown are assumed to be zero. As elucidated in Chen2008 [START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF], since the reference domain Laplacian model in multiple dimensions is not derived from an error equidistribution principle, its numerical solution until convergence is not strictly required to reach satisfactory mesh adaptation and, in practice, a number of Jacobi iterations in the order of O( 10) are generally su cient to reach the desired adaptation. The number of iterations n it will be reported for each case.

In section sec:cubeAn 4.1 adaptation is performed on a a steady Gaussian-like function, in order to test the convergence order on the interpolation error. In section sec:shockAn

adaptation is

performed on an unsteady analytical moving front passing over a sphere, in order to assess the capability of the model to preserve the validity of the mesh over intersecting curved boundaries throughout the time simulation.

Steady adaptation in a square and a cube sec:cubeAn

We consider the approximation of the function (↵, ↵ ) = (500, 0.1) in 3D. Also note that the a-posteriori limiter for the displacement is only applied in 3D, which is the case in which tangling is more often occurring.

⇢ = e ✓ 2 , = kx x xk 2 R 2 (20
On these meshes, we measure the L 2 error convergence of the P 1 interpolation ⇧⇢

||e|| L 2 = ✓Z ⌦ |⇢ ⇧⇢| 2 d⌦ ◆ 1 2 (21)
We plot the observed trends in figures fig: ErrorTrend2DConv 3a and fig: ErrorTrend3DConv 3b. It can be seen that in two dimensions it is easier to preserve, quite independently from the number of iterations n it performed, the second order convergence rate of the P 1 interpolation, with an error reduction for a given number of nodes shown in table

tab:Error2DConv
3, but a high number of iterations on a coarse mesh can actually increase the error.

In three dimensions, while the error on the adapted meshes is considerably lower (table

tab:Error3DConv 4
), the number of Jacobi iterations has to be increased to preserve the second order rate. Some adapted meshes obtained from the h = 0.1 and h = 0.05 initial meshes are visualized in figures possibly increasing the approximation error on the adapted mesh (as seen in table

tab:Error2DConv

3 for the 2D case for the coarsest meshes). The second consequence is that more iterations are needed on fine meshes to preserve the second order rate, as shown in figure

fig: ErrorTrend3DConv 3b. In three dimensions, the a-posteriori limiter also contributes to this e↵ect by constraining the allowed displacement of each node inside its ball at each iteration.

These e↵ects can be appreciated by observing the trend for the tetrahedron quality

Q = ⇣ P 6 j=1 l 2 j ⌘ 3/2 ↵|⌦ K | ( 22
)
where l j is the length of each edge of the element, |⌦ K | its volume, and ↵ the normalization factor to get Q = 1 on a regular tetrahedron with unit edges. Since r-adaptation inevitably introduces some anisotropy which is not taken into account in our quality measure, we expect the quality to be somewhat degraded in the adapted regions. Any-

[k] = (1 E [k] /E [0]
) with respect to the nonadapted case.

tab:Error2DConv h E [0] E [10] r [10] E [150]
[k] = (1 E [k] /E [0]
) with respect to the nonadapted case.

tab:Error3DConv

localized in the domain, can be a sign that the mesh is stretched also in smooth solution regions, possibly worsening the error reduction performances. In figure 4g corresponds to a significantly degradation of the elements quality for the h = 0.1 mesh, expecially when increasing the number of iterations, with more than 24% of elements having Q < 0.2 for 150 iterations, much higher than for the h = 0.05 (less than 10%).

#elements (%)

Mesh h = 0.1 The algorithm was tested by adapting over a moving front defined as

1.0 > Q > 0.8 0.8 > Q > 0.6 0.6 > Q > 0.4 0.4 > Q > 0.2 0.2 > Q > 0.0 1
⇢ (X(x, t)) = 8 > < > : 1 i f X(x, t) < 0 0.5 cos (s⇡X(x, t) + 1) if X(x, t) 2 [0, ] 0 i f X(x, t) > (23)
with

X(x, t) = s(x x 0 + vt) (24)
and scaling s = 20, initial position x 0 = 0.7, speed v = 0.2, front thickness = 0.005.

Unsteady mesh adaptation is performed on this analytical solution every t = 0.25.

The setup is shown in figures 6b. The domain is a quarter cylinder of radius 1.5 along the x-axis with x 2 [ 1.5, 1.5], surrounding a quarter sphere centered at the origin with radius 0.5. This case is designed to test as many geometrical sources of mesh tangling as possible before the application to fluid flow simulations, as it contains at the same time curved surfaces, ridges (the intersection of the sphere with each symmetry planes) and corners (the intersections of the sphere with both the symmetry planes), and a sharp solution moving over the geometry. Adaptation is performed with (↵, ↵ ) = (40, 0.1), with 30 Jacobi iterations, on an uniform mesh with edge size h = 0.05. The number of nodes and elements is reported in table tab:MeshDataFlow 5, as this is the same base mesh that will be used for the shock-sphere interaction simulations in the next section.

The obtained meshes are shown in figure Remarks on mesh folding and the purpose of the a-posteriori limiter. As discussed in section sec:Intro 1, there is no analytical proof for the validity of the meshes produced by our model neither in the continuum nor in the discrete setting. Examples of folded meshes have indeed already been reported in the literature for several other methods Dvinsky1991,Knupp1993 [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF][START_REF] Knupp | Fundamentals of grid generation, The Fundamentals of Grid Generation[END_REF].

Mesh folding has not been reported for the variable-di↵usion Laplacian in the reference domain in two dimensions

Ceniceros2001,Tang2003,Chen2008 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF][START_REF] Tang | Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws[END_REF][START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF], but in Ceniceros2001 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF] the authors themselves remark that there is no theoretical reason against its occurrence. In three dimensions, we have found that it is quite frequent to produce folded elements for too strong adaptation parameters or on concave boundaries when the limiter presented in the previous section is not applied. An example of the first situation is given in figure 7b, where two points on the surface are blocked and cannot move without folding the adjacent elements (the volume limiter is not applied, but displacement on the surface is limited on the surface ball in order to allow the projection on Bézier patches), and one element near the lower circle is folded.

In the numerical simulations presented in the next section, all of which have concave boundaries, tangling was observed whenever a shock wave hit or developed on the front of the object, without limiter. Since this happened in the first instants of the simulations, we have found that the straightforward three-dimensional extension of the original variable-di↵usion Laplacian method in the reference domain

Ceniceros2001

[10] would simply be unpractical on those cases without an additional limiting or correction step to avoid mesh folding.

Adaptation for unsteady compressible flows sec:FlowRes

We consider the simulation of unsteady inviscid compressible flows in a time dependent frame of reference. In particular, we couple the Fmg library we developed to the

Flowmesh solver

Guardone2011,Isola2015,Re2017 [START_REF] Guardone | Arbitrary lagrangian eulerian formulation for twodimensional flows using dynamic meshes with edge swapping[END_REF][START_REF] Isola | Finite-volume solution of two-dimensional compressible flows over dynamic adaptive grids[END_REF][START_REF] Re | An interpolation-free ale scheme for unsteady inviscid flows computations with large boundary displacements over three-dimensional adaptive grids[END_REF], based on a node-centered second order, total variationdiminishing finite volume scheme for the Euler equations, written in an Arbitrary-Lagrangian-Eulerian (ALE) form

Donea2004 [13] d dt Z ⌦(t) u d⌦ + I @⌦(t) n • (F(u) vu) d = 0 (25) eq:solvALE
where u is the array of the conservative solution, F(u) its flux, ⇢ denotes the mass density, ⇢U the momentum, and ⇢e t the total energy density. The moving domain velocity is represented by the vector field

v u = 0 @ ⇢ ⇢U ⇢e t 1 A T , F(u) = 0 @ ⇢U ⇢U ⌦ U + P I ⇢e t U + P U 1 A T (26) fig:FrontInputVol

fig:FrontTopCut

(f) Output mesh volumic cut at t = 3.5.

fig:FrontLeaveVol

(g) Output mesh boundary at t = 6.0.

fig:FrontLeaveCut

(h) Output mesh volumic cut at t = 6.0.

fig:Tangling

The pressure P is computed using the ideal gas equation of state for ideal gases

P = ( 1) ✓ ⇢e t 1 2 ⇢|U| 2 ◆
Within the code, a local conservative solution transfer procedure at each time step is guaranteed by the ALE formulation.

Unsteady mesh adaptation is performed according to the scheme shown in section eq:r:adaptDyn 3.5.

At each time step, the flow solution is predicted on the previous computational mesh, then the computational mesh is adapted, and finally the flow solution is recomputed on the adapted mesh. To this end, Flowmesh makes use of a conservative ALE-remap exactly matching the volumes swept by cell faces during mesh displacements and nodal volumes, and automatically fulfilling a Discrete Geometric Conservation Law (DGCL) GuillardFarhat2000,Farhat2001,Etienne2009 [START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF][START_REF] Farhat | The discrete geometric conservation law and the nonlinear stability of ale schemes for the solution of flow problems on moving grids[END_REF][START_REF] Étienne | Perspective on the geometric conservation law and finite element methods for ale simulations of incompressible flow[END_REF]. The code also includes the support of topological mesh modifications like edge split, edge collapse, barycentric node insertion, and Delauney node insertion, not used in this work.

To apply mesh adaptation at each time step, a low order computation of the solution at the next time step on the current mesh is used to provide a monitor function to the mesh PDEs.

Case 1: two-dimensional forward facing step

As a preliminary validation, we reproduce the results shown for the same method without a posteriori relaxation in Chen2008 [START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF] for the two-dimensional forward facing step fig: 3dcases we had to produce a refined mesh that is more than ten times bigger (in terms of nodes and elements) than the coarse one. Note that the diagonal cut is possibly the most demanding plane on which results can be compared, as the Laplacian model is uncoupled in multiple space directions, thus it tends to provide better results on cartesian planes, as shown in section tab:MeshTimeFlow 6. The benefits in terms of computational times in three dimensions are greater than in two dimensions. Anyway, while in two dimension we observed that mesh tangling was a rare occurrence with our Laplacian model, in three dimensions it was impossible to continue the time simulation without the a-posteriori limiter after the first few time steps, due to the strong deformation that quickly led to tangled elements at the step front and around its corners, but also at the shock reflection lines.

Case 3: shock-sphere interaction

In order to test the capabilities of the method to handle simultaneously shock waves and curved boundary, we choose to simulate the interaction of a traveling shock wave on a sphere. Some configurations for the di↵raction of shock waves over cylindrical and spherical obstacles have been studied experimentally for example in Bryson1961,Tanno2003 [START_REF] Bryson | Di↵raction of strong shocks by cones, cylinders, and spheres[END_REF][START_REF] Tanno | Interaction of a shock with a sphere suspended in a vertical shock tube[END_REF]. An early application of unstructured mesh adaptation to two-dimensional shock-cylinder simulations can be found in Drikakis1997 [START_REF] Drikakis | Computation of non-stationary shockwave/cylinder interaction using adaptive-grid methods[END_REF], while structured grid adaptation on axisymmetric shock-sphere simulations can be found in Sun2005 [START_REF] Sun | Unsteady drag on a sphere by shock wave loading[END_REF].

The simulation is limited to a quarter of a cylindrical domain (as for the analytically moving shock of the previous section, see figure

6.

In this case too it was impossible to complete the simulation over valid meshes without the action of the a-posteriori limiter near the corners and the curved surface.

Conclusions sec:Conclusions

The proposed algorithm for dynamic r-adaptation extends to three dimensions the method first proposed in Ceniceros2001,Chen2008,Arpaia2020 [START_REF] Ceniceros | An e cient dynamically adaptive mesh for potentially singular solutions[END_REF][START_REF] Chen | Second-order accurate godunov scheme for multicomponent flows on moving triangular meshes[END_REF][START_REF]Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF] for two-dimensional flows. An iterative solver based on diagonal Jacobi iterations for the discretized mesh PDEs with natural boundary conditions allows a cheap, uncoupled solution in each space direction. A novel a-posteriori relaxation scheme allows to prevent mesh tangling through the construction of a sequence of valid meshes also over curved boundary surfaces and corners, which is the main concern of r-adaptation methods in multiple dimensions, and it is interleaved with The reference domain formulation for mesh movement produces su ciently adapted meshes in as few as ten Jacobi iterations per time step during an unsteady flow simulation. While the a-posteriori relaxation algorithm is akin to a forward substitution algorithm, and thus formally dependent from the node ordering, this doesn't appear to spoil the adaptation pattern in any of our tests. We show the successful generation of valid adapted mesh on three-dimensional cases with moving shock waves. While the computational time overhead with respect to the original unadapted mesh is nonnegligeable, it is more than acceptable when compared to the simulation times needed to achieve the same accuracy on discontinuous flow features on uniformely refined mesh.

The attractiveness of the method rests in fact in its applicability on moving shocks, where an o↵-line mesh refinement approach would require to refine the mesh in most of the computational domain, and its easy coupling with ALE solvers, enabling solution conservation on the adapted meshes.

Limitations of this r-adaptation method are the same of the original two-dimensional formulation, namely the Laplacian models excessively pulls nodes towards non-convex boundaries and the displacement uncoupling in the multiple space directions can create sensible adaptation patterns for excessively strong adaptation parameters. Also, the choice of the parameters of the monitor function appear to be application dependent.

possibly leading to excessive mesh stretching for same values of the parameters. In these extreme situations, the e↵ect of the novel a-posteriori relaxation scheme allows nonetheless to recover a valid mesh by blocking mesh displacement in critical zones, allowing to continue the mesh movement at successive time steps as the flow features evolve away from the blocked mesh elements. We would like to remark that our limiting procedure is targeted at preserving mesh validity throughout the adaptation procedure.

This means that a nodal displacement can be blocked if the volume of an adjacent element falls below an user-defined threshold, but the mesh remains valid. Thus, the vertex positions of the blocked elements can be relaxed either by a subsequent application of r-adaptation at the next time step as the monitor function moves (as it is often the case in the simulation of traveling waves), either by the application of standard smoothing algorithms, which are fundamentally simpler than untangling methods.

While a linear finite element approximation is su cient to model the nodal degrees of freedom of straight-sided meshes, generalizations of the a-posteriori limiting method to curved meshes can be envisaged by increasing the degree of the finite element basis. This would require the formulation of a volume positivity predicate for the curved tetrahedron, which is outside the scope of this work.

The developments shown in this work have been primarily motivated by their application to capturing moving shock waves in inviscid compressible flows. However, r-adaptation can have interesting applications also to viscous flows, as preliminary work nouveau:tel-01500093 [START_REF] Nouveau | Adaptive residual based schemes for solving the penalized Navier Stokes equations with moving bodies : application to ice shedding trajectories[END_REF] shows that vorticity-based sensors allow capturing vortical structures in separated flows. This, together with the critical treatment of boundary layers in turbulent flow computations, will be investigated in the future. Future research lines also include the parallelization of the current method, for which no specific problems are envisaged, and the study of r-adaptation as a tool to complement h-adaptation in time-dependent simulations to somewhat reduce the overhead of the adaptation strategy.

  our current implementation, the re-evaluation is performed by linearly interpolating the scalar field f at the current nodal positions x [k] i through a standard search algorithm based on barycentric coordinates.

, 1 (

 1 fig:relax 1

Figure 1 :

 1 Figure 1: Two-dimensional illustration of the nodal displacement limiting. fig:relaxStart 1a: Proposed displacement for node i right after the Jacobi iteration k + 1, that would produce inverted elements (in red). fig:relaxCorrect 1b: Relaxed displacement producing valid elements (in green). fig:relaxOk 1c: Updated configuration. fig:relax

tion sec:NumModelBCs 2 . 1 ,

 21 slip boundary conditions bring a position constraint expressing the belonging of the node to the surface, and a null normal stress condition on the local tangent plane.The latter being already fulfilled by the Neumann conditions naturally imposed on the weak formulation, only the former is of our interest here. Although node belonging to an intersection curve can be formally expressed as the belonging to the two surfaces sharing the curve, this is not practical from an implementation point of view. Slip boundary conditions need thus to be specialized to the chosen geometry approximation and to distinguish among regular curved surfaces, ridges, and corners.In the following, the boundary treatment is detailed for the supported geometrical features: manifold surfaces, ridges (i.e. intersections of two manifold surfaces) and corners (intersections of two or more ridges). Note that di↵erent geometrical representations involving other local or global manifold parametrizations can be easily embedded in the algorithm.

Figure 2 :

 2 Figure 2: Illustration of the slip boundary projection procedure.

  fig:mesh-proj 2-(a)

) with ✓ = 40 , 2 .

 402 fig:solutionCube1 4a and fig:solutionCube2 4b. In both cases we consider a series of simplicial meshes with a uniform mesh size distribution, and di↵erent average edge size h, whose details are shown in tables tab:MeshData2DConv 1 and tab:MeshData3DConv 2. The above function is chosen in order to test capability of the models to adapt on a circle represented by a smooth solution field, before their application to solutions with sharp/discontinuous features. The mesh PDE parameters are set to (⌧, ⌧ ) = (5000, 1.0) in 2D, and to

  fig:cubesCut 4 to help understand these two phenomena. Taking as example the three-dimensional case, as the initial mesh is refined from h = 0.1 to h = 0.05 in figure fig:cubesCut 4, it can be appreciated that the displacement produced by the same number of iterations and the same adaptation parameters is smaller. This has two consequences. The first consequence is that a high number of iterations on coarse meshes can excessively stretch the mesh elements (as shown in figure fig:cubesCut010100 4g) in an orthogonal pattern, due to the uncoupling of the Laplacian model in the coordinate directions,

  fig:quality 5 we plot the evolution of the histograms of the elements quality with the number of iterations for the h = 0.1 and h = 0.05 meshes. The excessive stretch observed in figure fig:cubesCut010100

2 ⇠ h trend ⇠ h 2

 22 fig:ErrorTrend2DConv

2 ⇠ h trend ⇠ h 2

 22 fig:ErrorTrend3DConv

Figure 3 :

 3 Figure 3: Interpolation error convergence with mesh adaptation for the square and cube analytical test cases. fig:convergence

Figure 4 :

 4 Figure 4: Monitor function (top row) and volumic cuts in the adapted meshes (second to last row) for the cube test case, for di↵erent number of iterations, on the h = 0.1 and h = 0.05 initial meshes.fig:cubesCut

  fig:FrontInputVol 6a and fig:FrontInputCut

  fig:MovingFront 6 showing in particularly that the method is able to preserve a valid mesh both when the front is passing over the surface of the sphere (figures fig:FrontTopVol 6e, fig:FrontTopCut 6f) and most importantly when it hits and leaves the sphere

  fig:TanglingCube 7a, where an inverted element is produced just outside of the most refined region. An example of tangling on a concave boundary is given in figure fig:TanglingSphere

( a )

 a fig:FrontInputCut

Figure 6 :

 6 Figure 6: Moving front test case, meshes from t = 0.0 to t = 6.0.

Figure 7 :

 7 Figure 7: Examples of folded meshes with no limiter applied.

10 .[

 10 fig:contour_t1008

5. 2 .

 2 fig:gridinitstep 12a). Adaptation is performed on the mass density (figures fig:mass_step70??

Figure 8 :

 8 Figure 8: Two-dimensional forward facing step mass density contour lines at t = 0.5 and t = 1.0.

Figure 9 :

 9 Figure 9: Two-dimensional forward facing step adapted meshes at t = 0.5 and t = 1.0.

Figure 10 :

 10 Figure 10: Two-dimensional forward facing step mass density contour lines at t = 1.5 and t = 2.0.

Figure 11 :

 11 Figure 11: Two-dimensional forward facing step adapted meshes at t = 1.5 and t = 2.0.

Figure 12 :

 12 Figure 12: Initial meshes, adapted meshes and solution for the three-dimensional forward facing step and shock-sphere interaction cases.

  fig:meshstep114 and fig:meshstep216. Computational

  fig:gridinitsphere 12b). We choose a planar shock moving at M s = 1.5. Adaptation is performed on the mass density (figures fig:mass_sphere90 ?? and fig:mesh_sphere90 ??), with (↵, ↵ ) = (40, 0.1). Again, the aim is to compare the results obtained with mesh adaptation on a base mesh with edge size h = 0.05 with those obtained on a uniformly refined mesh with edge size h = 0.02. Mesh data are shown in table tab:MeshDataFlow 5. Contour lines

  (a) Nonadapted h = 0.04 mesh at t = 0.5. (b) Nonadapted h = 0.04 mesh at t = 1.0. (c) Adapted h = 0.04 mesh at t = 0.5. (d) Adapted h = 0.04 mesh at t = 1.0. (e) Fine (nonadapted) h = 0.015 mesh at t = 0.5. (f) Fine (nonadapted) h = 0.015 mesh at t = 1.0.

Figure 13 :

 13 Figure 13: Three-dimensional forward facing step mass density contour lines at t = 0.5 and t = 1.0.

Figure 14 :

 14 Figure 14: Three-dimensional forward facing step adapted meshes at t = 0.5 and t = 1.0.

Figure 15 :

 15 Figure 15: Three-dimensional forward facing step mass density contour lines at t = 1.5 and t = 2.0.

Figure 16 :

 16 Figure 16: Three-dimensional forward facing step adapted meshes at t = 1.5 and t = 2.0.

Figure 17 :

 17 Figure 17: Shock-sphere interaction mass density contour lines at t = 0.5 and t = 1.0.

Figure 18 :

 18 Figure 18: Shock-sphere interaction adapted meshes at t = 0.5 and t = 1.0.

Figure 19 :

 19 Figure 19: Shock-sphere interaction mass density contour lines at t = 1.5 and t = 2.0.

Figure 20 :

 20 Figure 20: Shock-sphere interaction adapted meshes at t = 1.5 and t = 2.0.

Table 1 :

 1 Mesh data for the 2D square convergence analysis.

	way, a too high percentage of bad quality elements, when sharp solution fronts are quite

Table 2 :

 2 Mesh data for the 3D cube convergence analysis.

	tab:MeshData3DConv
	h 0.15	E [0] 1.525974e-01 5.693340e-02 62.6905 % 1.194055e-01 21.751 % E [10] r [10] r [150] E [150]
	0.1	7.379553e-02 3.313339e-02 55.1011 % 4.372964e-02 40.742 %
	0.075 4.288518e-02 2.097308e-02 51.0948 % 1.991813e-02 53.555 %
	0.05	1.958636e-02 1.243068e-02 36.5340 % 1.090836e-02 44.306 %
	0.025 4.974168e-03 3.788809e-03 23.8303 % 3.614722e-03 27.330 %
	0.0125 1.258864e-03 1.142442e-03 9.2482 % 9.152757e-04 27.294 %

Table 3 :

 3 Interpolation errors E[k] = ||e[k] || L 2 for the 2D square convergence analysis, for 10 and 150 iterations, and reduction r

Table 4 :

 4 Interpolation errors E[k] = ||e[k] || L 2 for the 3D cube convergence analysis, for 10 and 150 iterations, and reduction r

  Evolution of mesh elements quality Q with the number of iterations nit for the 3D cube test, for the h = 0.1 and h = 0.05 meshes.

					Mesh h = 0.05
		100	1.21 4.43	1.35 5.41	1.53 6.54	1.78 7.71	2.08 8.88	5.45	9.78
	#elements (%)	40 60 80	13.61 40.3	15.9 40.1	17.61 19.09 20.57 39.76 39.44 38.89	14.18 24.47 35.43	15.53 25.31 33.64
		20	40.45 37.24 34.55 31.97 29.58
									20.48 15.74
			0	2 0	3 0	4 0	5 0	1 0 0	1 5 0
						n it
	Figure 5:

fig:quality 4.2. Moving front passing over a spherical boundary sec:shockAn

Table 5 :

 5 Number of nodes and elements for the simplicial meshes employed for the unsteady compressible flow cases.

	Emery1968,Leer1979,Woodward [17,

Table 6 :

 6 Computational times comparison. The overhead due to solution prediction and adaptation is important, but negligible if compared with an uniform refinement strategy.

Acknowledgements

We would like to recognize the fundamental contribution our colleague and friend Cécile Dobrzynski brought to this work.


Moya Jones

Framing Welsh identity

Keywords:

Welsh Studies as a cross-disciplinary field is growing both in Wales and beyond. Historians, sociologists, political scientists and others are increasingly collaborating in their study of the evolution of Welsh identity. This concept which was for so long monopolised and marked by a strong reference to the Welsh language and the importance of having an ethnic Welsh identity is now giving way to a more inclusive notion of what it means to be Welsh. Welsh identity has become a more open and inclusive notion since the institutional reforms brought about by devolution have anchored themselves in daily realities and subsequently in conceptual approaches to the problematizing of Welsh identity.

Résumé

Les études galloises constituent un nouveau champ inter-disciplinaire où historiens, sociologues, politologues etc., au pays de Galles même et au-delà de ses frontières, se rencontrent dans leur recherche sur l'identité galloise contemporaine. Les réformes postdévolutionnaires ont induit des changements dans les comportements et les mentalités des Gallois, entraînant une évolution dans la définition de l'identité galloise. Cette nouvelle identité n'est plus définie par les critères ethno-linguistiques d'antan mais se situe plutôt dans une nouvelle aire où les notions de citoyenneté et d'inclusion priment.

America) while within Wales researchers and intellectuals are beginning to leave their own well-defined fields of study (geography, history, literature, sociology, politics, gender studies …) to engage in more cross-disciplinary projects and programmes. These collaborative efforts are concentrated on investigations into Wales and Welsh identity/identities. Framing Welsh identity is thus a question which is exercising many people at present as the Welsh Assembly, by its very existence, is creating and defining its own legitimacy. One of the basic interrogations underpinning this surge of activity is the question of whether today one refers to « Welsh people » or the « people of Wales», for the tricky problem of identity based on language skills is still very much to the fore. These aspects do not arise in the case of Scotland and Scottish identity because Scotland's past as a unified kingdom, the maintained presence of its institutions, the absence of a major language problem and the strong nationalist demand for devolution make it a clearly different case to Wales. Who is Welsh? The two intertwined, interpenetrative elements of territory and language have fashioned Welsh identity, as has the country's resulting history. Difficult terrain fashioned communities and mentalities while hindering unity. The nature of the land also determined industrial and agricultural differences and subsequent political divisions; thus the pattern and nature of nationalism and politics in Wales have hinged on factors that have been in place for centuries. Wales covers about 21,000 sq km which represents about a tenth of the British land mass although it is home to only about 2.9 million people, roughly 5% of the total UK population. The population is concentrated along the flat coastlands especially in the south around the cities of Cardiff, the capital, Swansea and Newport, and in the narrow valleys that run down to the sea. These areas were long the site of industrial activity such as metallurgy and mining. The central high plateau, dominated by rural activities, has always made communication and transport difficult within Wales. In the centre and the north-east the low population density has played a crucial role in the survival of the language. Similarly, communication difficulties led to comparatively isolated communities being self-sufficient and even inward looking, receptive to the ideas of self help and independence fostered by Nonconformist religion, very close knit, distinguishing themselves clearly from other social units like other villages or religious groups. This tendency to identify with smaller rather than larger entities is still alive today, reflected internally in identification with a neighbourhood rather than with a town, and in Wales's outward policies, reflected in European regional affiliations. Geography also explains in part why Wales never experienced the political unity that the kingdom of Scotland once enjoyed. Gwynfor Evans, the leading figure of modern Welsh nationalism, blamed geography and the resulting historical and political scattering of the Welsh across their country for a lack of political and national cohesion : « Ironically, one of the reasons Wales has found it difficult to obtain decentralisation is that the Welsh themselves are too decentralist to unite » (Evans 2003 : 215) However, we could take the view that it is History that gives us the best view to framing Welsh identity for there has been a striking and sustained renaissance in the study of Welsh history over the past two generations. There was a dearth of research in this area until the 1960s when a number of authoritative historians such as Glanmor Williams, K O Morgan, and Gwyn Alf Williams at Swansea and Aberystwyth undertook a transformation of this field. These first historians were preoccupied with a need to highlight national unity and the wholeness of Wales, attempting to pinpoint « When was Wales? » [START_REF] Williams | When Was Wales? A History of the Welsh[END_REF]. Since then, institutional and social factors have prompted the rise of a group of historians who have mainly been concerned with social history and empirical, archival research, writing a history 'from below'. Some have written more from a committed, political viewpoint but what connects such eminent researchers as Chris Williams, Geraint Jenkins, Duncan Tanner, Bill Jones and John Davies is their view of themselves as « Remembrancers » (« passeurs de mémoire »). However, while inevitably the actions of forces both inside and outside Wales (globalization, mass media, movement of people, information technology …) are transforming Wales into a more mutable and fragmented community of peoples, Welsh history cannot be fully understood without reference to its relationship with England. The story of Wales is not so much the history of kings as the history of the people, « y gwerin », and of their political and cultural leaders. Since the first incursions of the English King Edward I in the late 13th century Wales has been subsumed into the history of England (the Encyclopaedia Britannica 1888 entry reads: « For ' Wales' see 'England' »), and Gwynfor Evans was adamant that ignorance of their own history was -and is -the reason why many Welsh people have been slow in coming to a sense of nationhood: « In Wales, the weakened sense of a collective past, and, therefore, of a collective future, is not only due to a lack of national institutions. » (Evans 1996 : 134). His guidance and leadership of Plaid Cymru, the Welsh nationalist party, were to be driven by a strong sense of the long term cultural unity which has held the country together longer and better than any political arrangement. The joined histories of England and Wales have been marked, since the Acts of Union in 1536 and 1542, by a symbiotic relationship that has at times been mutually beneficial but which many believe to have been mainly dominated by England to the detriment of Wales and its distinct identity. Establishing English law, excluding Welsh speakers from public office and enforcing the supremacy of Church of England were the first instances of this. A distinct Welsh identity continued to exist in matters of religion, with the growth of Puritanism in Wales, especially Methodism in the 18th century, and the Welsh language was often the vehicle of that faith. However, the Industrial Revolution of the 19 th century brought an influx of Anglophone immigrants mainly to the south of the country thus contributing to the decline of the language and to a certain dilution of identity. Politically, Wales had been dominated by the Liberal party and the Radical movement but allegiance to the Labour Party and class solidarity became the political characteristics of the country during most of the 20 th century. Economic depression in the 1930s and the post-war decline of industry (of which the 1984-1985 miners' strike is perhaps the most striking event) has transformed the country. Economic structures have changed as deindustrialisation has been accompanied by the arrival of light industries, service industries, the feminisation of the workforce and more prominent inward investment by British and overseas concerns. Also, religious observance has declined, as it has in the rest of Britain and this has been in parallel with the erosion of the Welsh language.

Gwnewch bopeth yn Gymraeg ! (do everything in Welsh).

Today approximately 600,000 people speak Welsh, which is about 20% of the population and the table below shows how the number of Welsh speakers has declined over the years. This decline is easily explained by industrialization, the proximity of England and indeed globalization, and it has accelerated as modern forms of communication, TV, radio, press, and the internet have expanded. Is language therefore an element to be taken into account when attempting to define Welsh identity? The early struggle for recognition of Wales as a separate country or region and the first claims for autonomous government, and indeed the founding of Plaid Cymru, the nationalist party of Wales (1925) were based on what could be termed as mythical values surrounding Welsh cultural heritage which was intimately linked to the language. Modern consciousness of the importance of the language as an identity marker and of the danger threatened by its decline was first aroused by the poet Saunders Lewis in a 1965 BBC radio broadcast -« Tynged yr Iaith » (The Fate of the Language) warning of the disappearance of the language, and of specific Welsh identity. A long and sometimes violent campaign by militants from Cymdeithas yr Iaith Cymraeg (The Welsh Language Society) eventually led to the passing of two laws (1967 Welsh Language Act and the 1993 Welsh Language Act) to protect and promote the language, giving it parity with English, and making the learning of Welsh obligatory in schools. The recent increase in the number of speakers is due to efforts by the Welsh Language Board in conjunction with the Welsh Assembly's language promotion programme, Iaith Pawb (Language for All). The defenders of the language promote it as giving «added value» to the country, to Welsh identity and argue that the language enjoys colossal good will and support. It is hoped that the new awareness of all things Welsh, the growth in the awareness of national identity that seems to have accompanied devolution will also promote Welsh cultural identity and the movement is inspired by the example of other countries like Finland which, in the wake of political autonomy, have come to enjoy a cultural and linguistic reinforcement. Their dream is of a truly bi-lingual Wales. This project is supported by the fact that the UK government has ratified the 1992 Council of Europe European Charter for Regional or Minority Languages.

The detractors believe that it is useless to promote the language as so few people really speak it. Often described as a « barbaric » language, and attacked throughout history, from the Acts of Union to the notorious « Blue Books » (the 1847 Education Report) the language is still a divisive element today among the Welsh themselves. Welsh speakers are viewed with hostility by Anglophone Welsh for they view Welsh speakers as a self-appointed elite who sideline non Welsh speakers and who reserve promotions and favours within Wales and Welsh institutions for a select few. Nevertheless, despite the criticism of the way statistics are defined and gathered, we can confidently say that Welsh is the everyday language of about half a million people throughout Wales although its use can be tied clearly to certain areas. In 1985 Denis Balsom determined that Wales could be divided into three main language areas, according to self-identification (which has since been enhanced and modified by other researchers). Those zones were, firstly, « British Wales », which was the eastern half of Wales, the Cardiff area and Pembroke in the extreme south-west. In this zone people identified with Britain and spoke English and no Welsh. « Welsh Wales » incorporated the rest of the south-east and included the Valleys. These people identified themselves with Wales but generally did not speak Welsh. Finally, « Y Fro Gymraeg » corresponded mainly to the western half of Wales and it is the area where the population identify with Wales and also speak Welsh. The last two areas have shared similar voting patterns in the 1997 referendum and in Assembly elections (Williams 2003 : 10). The language is today gaining ground and more speakers in urban areas such as Cardiff and Swansea and in Monmouthshire and Glamorgan. It is most under threat in the north and the west due to tourists, English retirees and second home buyers and to the cities of Liverpool and Manchester social services' « dumping » of problem families in the struggling guesthouses on north coast resorts like Rhyl and Llandudno.

Politics in Wales.

The language issue is inextricably linked to politics and in particular to the story of Plaid Cymru. The party was however long associated with a romantic, idealized view of an agrarian past, with a subsequently strong ethnic approach to the question of Welshness. In fact the party's problem in capturing votes was that it was (and to some extent still is) seen as « too Welsh » and therefore off-putting to some voters. Long associated with Gwynfor Evans, it sprang up in the line of a certain renaissance and a revival of Welsh sentiment, the view of Wales as a nation. Its demands for autonomy eventually shook the statewide Conservative and Labour parties, and following Gwynfor Evans's historic election as the nationalist MP for Carmarthen in 1966 the first moves were made to mollify the nationalists. We can say briefly that the first ideas for devolution were made by the Labour government in the late 60s in a context of constitutional reform but also to save Labour seats in Wales (and Scotland). There is some truth then in saying that, politically, the Labour party has also been instrumental in shaping Welsh identity. This happened initially through the support that Labour received in Wales when the Liberal party lost ground to it and when Labour came to represent class solidarity above all, which persisted throughout the latter half of the 20 th century. However, there may also be some truth in saying that since 1997 New Labour has been instrumental in helping to carve out a new identity for the Welsh. Comparing the results of the two referendums held in Wales helps to remind us of how little enthusiasm there was for the proposals set out in the 1978 Wales Act : Some of the reasons for this defeat in Wales include the fact that the amendment requiring 40% of the electorate to approve of the contents of the Wales Act and the fact that this legislation was framed beforehand. The economic and social climate of the time (the 1978-79 «winter of discontent») created a climate of mistrust in the government. Issues had become very polarized and state-wide there was dominating interest for Scotland. Fears of what devolution would entail included the break up of Britain and while Anglophones feared domination by Welsh speakers, those in the Welsh-speaking north were disgruntled by the idea that politics would be dominated by Anglophones in Cardiff. In any case, the Welsh context was diluted by British party politics, especially by internal Labour wrangling. Prominent anti-devolutionists included the future party leader, Neil Kinnock. However, eighteen years later British politics, Welsh politics and the population of Wales had evolved and changed.

The "quiet earthquake".

What had happened between the two referendums to explain the difference in the outcome? The reasons are various and complex but are deemed to have been influenced by the earlier, favourable vote in Scotland on 11 September 1997. While the devolution measures were virtually identical to those of 1979, the positive outcome has often been attributed to eighteen years of Conservative government in Westminster while support for the Tories dwindled in Wales. However, this period witnessed the rise in the number and prominence of institutions put in place by both Labour and Conservative governments which were allowing the development of a civic national identity in Wales which had created an increasingly strong «proto-State». Some of these institutions are quite venerable, such as the National Library and the University (1907). Cardiff was officially designated as the capital city in 1955, the Welsh Office was created in 1964 and the Welsh Development Agency in 1976. Meanwhile other bodies sprang up such as the Welsh Arts Council, the Wales Tourist Board, the Sports Council for Wales, the Welsh TUC and the Welsh CBI, the Welsh Football Association, the Welsh Rugby Union, the Welsh National Opera, etc., etc. Although Wales never knew a constitutional convention like Scotland these institutional creations have helped to create and develop a sense of nationhood. However, and perhaps just as importantly, the 1997 referendum results also point to the fact that a generational shift had taken place. A new population voted, reflecting demographic, economic and even psychological changes. Those under the age of 45 voted «yes» by a margin of 3 : 2 (Osmond 2007 : 20). Although there were 87,000 fewer voters in the second referendum there was a massive shift in the number of « yes» votes, up by 316,365 on the 1979 result while the «no» votes dropped by 404,632. Was this in fact a vote of confidence in Wales and Welshness and an affirmation of a new Welsh identity? John Osmond thinks so especially as regular polls have revealed that since 1997 the constitutional preference of a growing number of Welsh people is for a Welsh Parliament (1997 :19.6%;1999: 29.9%;2001 : 38.8%;2003 : 37.8%;[START_REF] Evans | We are Welsh, but…' -national identification on the margin[END_REF]: 42.1%. Osmond, 2007 : 21) : 21). The creation of the aforementioned institutions may have gone some way towards helping to shape a distinct Welsh identity within the British framework but an important characteristic of many of them was that they were « quangos » or quasi-autonomous non-governmental organizations, which created a series of problems because they were run by appointees or co-opted members in unelected bodies that were unaccountable to the electorate. The resulting « democracy deficit » contributed to the specific political atmosphere of the long period of Conservative government (1979 -1997). « The case for a measure of democratic control over Welsh life by the people of Wales was energized not retarded by the political, industrial and economic traumas of the Thatcherite years » (Williams 2000 : 241) Margaret Thatcher's affirmation that « there is no such thing as society » (31 October 1987) infuriated the communitarian values of the Welsh, forced women in particular out of the shadows to tell their stories and prompted the Welsh historians referred to earlier, and others, into more fruitful approaches to their discipline and research. Thus history, like politics, has been caught up in issues of religion, language and nationality in Wales and in the current debate over identity or rather identities. Perhaps the plural should be preferred because many feel and know today that there is no longer a homogenous Welsh identity. Since the end of the 1980s, in Wales there has arisen a kind of consensus, a move away from previous, narrow definitions of Welsh identity based on ethnic criteria of language and culture, which had maintained divisions within Welsh society, to a new kind of collective identity: The Welsh are in the process of being defined not in terms of shared occupational experience or common religious experience or the survival of an ancient European language or for contributing to the Welsh radical tradition, but rather by reference to the institutions that they inhabit, influence or react to. This new identity may lack the ethical and political imperatives that characterized Welsh life for two centuries but it increasingly appears to be the only identity available. (Williams 2000 : 240). It was becoming clear that many of the characteristics that marked out the people of Wales were competing and conflictive and they were gradually being unmasked as anachronistic and misleading. Thus, if a new identity was to be generated that could accommodate all those who had made Wales their home, all those who, whatever their birthplace, ethnic background, linguistic ability or even self definition, constituted the « Welsh people » it would and will have to be pluralist, inclusive and constructed on a basis of citizenship rather than on the possession of certain « essential » attributes, attitudes or beliefs. Wales needs to establish clearly a new civic identity as opposed to the old ethnic identity. Similarly Plaid Cymru had to rethink its image and its policies because the people of Wales could not recognize themselves in the Party's original, narrow, language based definitions. This new conception of identity was what underpinned Ron Davies' view of devolution. As Labour's first post-devolution Secretary of State for Wales he was known as «the architect of devolution» and had famously stated that « devolution is a process not an event » [START_REF] Davies | Devolution : a Process not an Event[END_REF]. Ron Davies further articulated the concept that Welsh identity was mutating and moving forward, emphasizing the fact that it had to be an inclusive identity, for all the people who lived in Wales, whatever their origin or background (Davies 1999 : 6-7). This very recent, shifting nature of Welsh civic and cultural identities goes some way to explaining changes in popular and political attitudes to devolution from 1979-97 and why and how language ceased to be such a decisive issue. But even at the present moment the debate is still moving and shifting in focus and not just in the world of politics or in academia. In the world of culture too Wales is re-shaping itself, redefining itself, re-inventing itself, partly in reaction to the new politics. There is a wide cultural renaissance within post devolution Wales that is epitomized in Peter Lord's Imaging the Nation (2000). In the three volumes of this work he set out to establish and recover the visual heritage of the country as a narrative of Wales. For Peter Lord, and others, the notion of cultural regeneration accompanies that of political rebirth : It has often been alleged that a national consciousness heavily conditioned by the needs of differentiation from a dominant neighbour is a characteristic Welsh weakness. As a result of the political and economic decline of that neighbour, a complementary growth in our own self-confidence and a wider change in perceptions of nationality, it is the hope of many at the beginning of the twenty-first century that this essentially colonized state of mind may at last be transcended. (Lord 2000 : 9) This observation is particularly interesting because of the reference to « colonialism » because academic discussion has recently embarked on a debate, on Welsh identity and problematizing Wales as to whether or not Wales can be qualified as a «post-colonial» country or not. Can it be qualified as a colony of England and is it now in a period after colonialism, witnessing the withdrawal of imperial forces? Or is Wales a case for «postcolonial» studies, where disparate forms and representations, reading practices and values circulate across the barrier between colonial rule and national independence? [START_REF] Aaron | Postcolonial Wales[END_REF]. Thus politicians, historians and cultural commentators are coming together in agreement that a new Welsh identity should accommodate all those who made Wales their home whatever their birthplace, ethnic origin, linguistic competence or even their self defined identification. The basis of citizenship should be pluralistic and inclusive and not depend on possession of essential (internal) attributes, attitudes or beliefs. This reflection is also accompanied by the realisation that the concept of the nation-state is now out of date and is no longer a viable goal for those who would wish to see a greater autonomy for Wales. Greater visibility within a European framework is deemed by many to be a more realistically achievable aim, given the fuzzy, unclear nature of Wales's borders -internal and external -and the implications of the very close relations with England, both the recognized and the unacknowledged.

A new identity.

Nevertheless, despite this progressive and dynamic reflection and investigation, the answer to the original question, « Who is Welsh? » still remains elusive. Those who have carried out enquiries into Welsh identity have found it shifting and difficult to pin down. It is less positively experienced than in Scotland. The concept of national identity was one examined by the Welsh Election Study Survey in 1979 which asked respondents to identify themselves as «British, Welsh, English or whatever» and they found that among those who voted «Yes» in the referendum 76% of them identified themselves as Welsh and 19% said they were British (Foulkes et al. 1983 : 210). «No» voters were 50% Welsh and 45% British. However, among those who declared they were normally Labour supporters and also Welsh identifiers only 33% voted in favour of the referendum proposal. This result shows that in 1979 class solidarity, expressed in the support for Labour, was stronger than a Welsh national identity. However, by the time the 1997 referendum was held more sophisticated analyses showed rather different results. Earlier in the May 1997 general election the people of Wales had elected candidates exclusively from the three parties committed to some measure of devolution and no Conservative MPs were elected to Westminster. It was discovered that most « Yes » votes came from the Western, Welsh-speaking areas and the Welsh-thinking former industrial Valleys. Also, young people had voted in favour of the Welsh Assembly, but it is not sure whether this indicated a growing real consciousness of Welshness or not. Interestingly, a large number of «Yes» votes were recorded in Monmouthshire, the county on the English border which has been attached to both England and Wales at different periods. Further research carried out in Wales 1997, 1999 and 2001, using the mobile « Moreno » scale of self-identification, revealed that a growing number of people identified themselves as Welsh, (though there was no drop in the number who also defined themselves as British).In 1997, 42.9% identified themselves as « Welsh » or « more Welsh than British » but by 2001 this figure had risen to 48.1%. [START_REF] Jones | National Identity in Post-Devolution Wales» in Devolution : What Difference Has It Made?[END_REF]. This seems to indicate a polarisation of national identity, which is stronger among younger generations. It also points to the fact that Welshness is becoming more attractive. However, research carried out in 1998, just before the first Welsh Assembly elections, also revealed some interesting results. Richard Haesly tried to discover how national identities and « imagined communities » [START_REF] Anderson | Imagined Communities. Reflections on the Origin and Spread of Nationalism[END_REF] are linked [START_REF] Haesly | Identifying Scotland and Wales : types of Scottish and Welsh identities[END_REF]. According to his findings, the Welsh imagined community appears to be fundamentally more contested than the more easily forged Scottish identity. He determined that the Welsh respondents were not nationalists, and not sure if Wales is a meaningful nation before concluding that Wales is seen as a divided nation, with a thin sense of identity and a sense that Welsh history is relatively unimportant. However, the interviewees did share a strong emotional bond with Wales that replaced rather than enhanced specific ideas about the subjective meaning of contemporary Welshness. This focussed on the importance of sport and the national teams, beating England and a general anti-Englishness. These people are sure they are Welsh, but not sure what precisely Wales is or who is Welsh. He discovered a pride in the culture, but which was (or is) rather low level. Further research carried out in Flintshire in north-east Wales revealed hierarchical notions of Welshness with respondents divulging perceived degrees of «real» Welshness, [START_REF] Evans | We are Welsh, but…' -national identification on the margin[END_REF]. Respondents self-identified as being «Welsh» but not as Welsh as someone from Aberystwyth (on the west coast) for example. The respondents seemed to acknowledge that their Welshness was an identity that was internalised but it was something that was constantly being negotiated as if living in a place on the margin of national space seems to shape conceptions of national identity in a modern national polity that is still in its early days. Contrast this with the general impression that throughout Wales there is nevertheless a very strong sense of « community » that binds the nation together (Le Disez 2005).

Conclusion.

Can Welsh identity be « framed » therefore? Can it in fact be pinned down or is it too diffuse and indeed multiple? Devolution and national government for Wales is helping the move towards a reinforced and inclusive sense of identity which is more specifically Welsh. Wales is shifting and changing and reinventing itself, almost despite the Welsh themselves sometimes for theexistence of a nationalistic national identity is not the same as it is in Scotland. The absence of such a type in Wales suggests that Scotland and Wales each represent a very different political environment with a differing national consciousness. If such a new, inclusive multi-faceted Welsh identity exists today then who can frame it? As Raymond Williams put it : « Who speaks for Wales? » [START_REF] Williams | Who Speaks for Wales?[END_REF]. Ten years after devolution and the creation of an all-Wales institution whose powers are slowly increasing, it has to be the National Assembly that assumes this task. Wales no longer possesses a homogenous identity and the old identity markers (language, religion…) are crumbling and disappearing. Welsh identity is becoming diverse, fragmented and multi-faceted : 2.1% of the population or 62,000 people are from ethnic minorities (2001 Census) and only the National Assembly can now hold the country together and represent Wales to and for the Welsh. In the face of globalisation and the internet, it is becoming increasingly harder to maintain cultural or national distinctiveness but Wales is shifting, changing and re-inventing itself. At a time when the Assembly is set to enjoy new legislative powers there are encouraging signs that the people of Wales are understanding the potential of the Assembly. Participation in the 2007 elections rose to 43.5 % from the low of 38% in 2003 (The Electoral Commission, 2007 www.electoralcommission.co.uk) and the voters elected a coalition government where Labour now share power with Plaid Cymru. Unprecedented and promising.

Table 1 .

 1 Numbers of Welsh speakers

	year	population	% Welsh speakers
	1891	1,519,000	54,4%
	1911	2,421,000	43,5%
	1931	2,593,000	36,8%
	1961	2,644,000	26%
	1991	2,891,000	18,7%
	2001	2,903,085	21%
	Source :Welsh Language Board, UK Office for National Statistics

Table 2 Wales

 2

	referendum results 1 March 1979	turnout 58.8%

Table 3 Wales

 3

	referendum results 18 September 1997	turnout 50.1%
	I agree that there should be a Welsh Assembly	559,419 votes	50.3% voters
	I do not agree that there should be a Welsh Assembly	552,698	49.7%
	Source : British election survey/ BBC Political research unit.