Keywords
Computer Science and Digital Science
 A6.2. Scientific computing, Numerical Analysis & Optimization
 A6.2.7. High performance computing
 A6.2.8. Computational geometry and meshes
 A6.5.1. Solid mechanics
 A6.5.2. Fluid mechanics
Other Research Topics and Application Domains
 B5.2.3. Aviation
 B5.2.4. Aerospace
 B9.5.1. Computer science
 B9.5.2. Mathematics
 B9.5.3. Physics
 B9.5.5. Mechanics
1 Team members, visitors, external collaborators
Research Scientists
 Frederic Alauzet [Team leader, Inria, Senior Researcher, HDR]
 PaulLouis George [Inria, Emeritus]
 Adrien Loseille [Inria, Researcher, HDR]
 David Marcum [Inria, International Chair, Advanced Research Position]
PhD Students
 Sofiane Benzait [CEA]
 Francesco Clerici [Inria]
 Lucien Rochery [Inria]
 Lucille Marie Tenkes [Inria]
Technical Staff
 Matthieu Maunoury [Inria, Engineer]
 Cosimo Tarsia Morisco [Inria, Engineer, from Mar 2021]
Administrative Assistant
 Maria Agustina Ronco [Inria]
External Collaborators
 Rémi Feuillet [SIEMENS INDUSTRY SOFTWARE, until Sep 2021]
 Loic Marechal [Dassault Systemes]
 Julien Vanharen [ONERA]
2 Overall objectives
Numerical simulation has been booming over the last thirty years, thanks to increasingly powerful numerical methods, computeraided design (CAD) and the mesh generation for complex 3D geometries, and the coming of supercomputers (HPC). The discipline is now mature and has become an integral part of design in science and engineering applications. This new status has lead scientists and engineers to consider numerical simulation of problems with ever increasing geometrical and physical complexities. A simple observation of this chart
shows: no mesh = no simulation along with "bad" mesh = wrong simulation. We have concluded that the mesh is at the core of the classical computational pipeline and a key component to significant improvements. Therefore, the requirements on meshing methods are an ever increasing need, with increased difficulty, to produce high quality meshes to enable reliable solution output predictions in an automated manner. These requirements on meshing or equivalent technologies cannot be removed and all approaches face similar issues.
In this context, Gamma team was created in 1996 and has focused on the development of robust automated mesh generation methods in 3D, which was clearly a bottleneck at that time when most of the numerical simulations were 2D. The team has been very successful in tetrahedral meshing with the wellknown software Ghs3d34, 35 which has been distributed worldwide so far and in hexahedral meshing with the software Hexotic55, 56 which was the first automated full hex mesher. The team has also worked on surface meshers with Yams25 and BLSurf15 and visualization with Medit. Before Medit, we were unable to visualize in real time 3D meshes !
In 2010, Gamma3 team has replaced Gamma with the choice to focus more on meshing for numerical simulations. The main goal was to emphasize and to strengthen the link between meshing technologies and numerical methods (flow or structure solvers). The metricbased anisotropic mesh adaptation strategy has been very successful with the development of many error estimates, the generation of highly anisotropic meshes, its application to compressible Euler and NavierStokes equations 8, and its extension to unsteady problems with moving geometries 11 leading to the development of several softwares Feflo.a/AMGLib, Wolf, Metrix, WolfInterpol. A significant accomplishment was the highfidelity prediction of the sonic boom emitted by supersonic aircraft 9. We were the first to compute a certified aircraft sonic boom propagation in the atmosphere, thanks to mesh adaptation. The team has started to work on parallelism with the development of the multithread library LPlib and the efficient management of memory using space filling curves, and the generation of large meshes (a billion of elements) 49. Theoretical work on highorder meshes has been also done 33.
Today, numerical simulation is an integral part of design in engineering applications with the main goal of reducing costs and speeding up the process of creating new design. Four main issues for industry are:
 Generation of a discrete surface mesh from a continuous CAD is the last nonautomated step of the design pipeline and, thus, the most human time consuming
 Highperformance computing (HPC) for all tools included in the design loop
 The cost in euros of a numerical simulation
 Certification of highfidelity numerical simulations by controlling errors and uncertainties.
Let us now discuss in more details each of these issues.
Generating a discrete surface mesh from a CAD geometry definition has been the numerical analysis Achille's heel for the last 30 years. Significant issues are far too common and range from persistent translation issues between systems that can produce ill defined geometry definitions to overwhelming complexity for full configurations with all components. A geometry definition that is ill defined often does not perfectly capture the geometry's features and leads to a bad mesh and a broken simulation. Unfortunately, CAD system design is essentially decoupled from the needs of numerical simulation and is largely driven by the those of manufacturing and other areas. As a result, this step of the numerical simulation pipeline is still labor intensive and the most time consuming. There is a need to develop alternative geometry processes and models that are more suitable for numerical simulations.
Companies working on hightech projects with high added value (Boeing, Safran, DassaultAviation, Ariane Group, ...) consider their design pipeline inside a HPC framework. Indeed, they are performing complex numerical simulations on complex geometries on a dailybasis, and they aim at using this in a shapeoptimization loop. Therefore, any tools added to their numerical platform should be HPC compliant. This means that all developments should consider hybrid parallelism, i.e., to be compatible with distributed memory architecture (MPI) and shared memory architecture (multithreaded), to achieve scalable parallelism.
One of the main goals of numerical simulation is to reduce the cost of creating new designs (e.g reduce the number of windtunnel and flight tests in the aircraft industry). The emergence of 3D printers is, in some cases, making tests easier to perform, faster and cheaper. It is thus mandatory to control the cost of the numerical simulations, in other word, it is important to use less resources to achieve the same accuracy. The cost takes into account the engineer time as well as the computing resources needed to perform the numerical simulation. The cost for one simulation can vary from 15 euros for simple models (1D2D), to 150 euros for Reynoldsaveraged NavierStokes (3D) stationary models, or up to 15 000 euros for unsteady models like LES or LatticeBoltzmann 1. It is important to know that a design loop is equivalent to performing between 100 and 1 000 numerical simulations. Consequently, the need for more efficient algorithms and processes is still a key factor.
Another crucial point is checking and certification of errors and uncertainties in highfidelity numerical simulations. These errors can come from several sources:

i)
modeling error (for example via turbulence models or initial conditions),

ii)
discretization error (due to the mesh),

iii)
geometry error (due to the representation of the design) and

iv)
implementation errors in the considered software.
The error assessment and mesh generation procedure employed in the aerospace industry for CFD simulations relies heavily on the experience of the CFD user. The inadequacy of this practice even for geometries frequently encountered in engineering practice has been highlighted in studies of the AIAA 2 CFD Drag Prediction Workshops 59 and HighLift Prediction Workshops 74, 73. These studies suggest that the range of scales present in the turbulent flow cannot be adequately resolved using meshes generated following what is considered best present practices. In this regard, anisotropic mesh adaptation is considered as the future, as stated in the NASA report "CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences" 76 and the study dedicated to mesh adaptation 65.
These preoccupations are the core of the Gamma project scientific program. To answer the first issue, Gamma will focus on designing and developing a geometry modeling framework specifically intended for mesh generation and numerical simulation purposes. This is a mandatory step for automated geometrymesh and mesh adaptation processes with an integrated geometry model. To answer the last three issues, the Gamma team will work on the development of a highorder meshadaptive solution platform compatible with HPC environment. To this end, Gamma will pursue its work on advanced mesh generation methods which should fulfill the following capabilities:

i)
geometric adaptive,

ii)
solution adaptive,

iii)
highorder,

iv)
multielements (structured or not), and

v)
using hybrid scalable parallelism.
Note that items $i)$ to $iv)$ are based on the wellposed metricbased theoretical framework. Moreover, Gamma will continue to work on robust flow solvers, solving the turbulent NavierStokes equations from second order using Finite Volume  Finite Element numerical scheme to higherorder using Flux Reconstruction (FR) method.
The combination of adaptation  highorder  multielements coupled with appropriate error estimates is for the team the way to go to reduce the cost of numerical simulations while ensuring highfidelity in a fully automated framework.
3 Research program
The main axes are:
 Geometric Modeling:
 Highfidelity discrete CAD kernel.
 Continuous parametric CAD kernel.
 Enhanced Generic Meshing Algorithm:
 Adaptation (extreme anisotropy, metricaligned, metricorthogonal).
 Highorder (tetrahedra, hexahedra, boundary layer, adapted).
 Larges meshes (tetrahedra, hexahedra, adapted).
 Moving mesh methods for moving geometries.
 Toward Certified Solutions to the NavierStokes Equations:
 Flow solver and adjoints (Finite Volumes, Finite Elements, Flux Reconstruction).
 Error estimates and correctors.
 Advanced Mesh and Solution Visualisation:
 Pixel exact rendering (HighOrder mesh, HighOrder solution).
 Preprocessing and postprocessing.
4 Application domains
Our research in mesh generation, mesh adaptation and certification of the Numerical Simulation Pipeline finds applications in several different domains such as aviation and aerospace but also all fields where computation and simulation are used: fluid mechanics, solid mechanics, solving wave equations (acoustic, electromagnetism...), energy or biomedical.
5 New software and platforms
5.1 New software
5.1.1 FEFLOAREMESH

Keywords:
Scientific calculation, Anisotropic, Mesh adaptation

Functional Description:
FEFLOAREMESH is intended to generate adapted 2D, surface and volume meshes by using a unique cavitybased operator. The metricaligned or metricorthogonal approach is used to generate high quality surface and volume meshes independently of the anisotropy involved.
 URL:

Contact:
Adrien Loseille

Participants:
Adrien Loseille, Frederic Alauzet, Rémi Feuillet, Lucien Rochery, Lucille Marie Tenkes
5.1.2 GHS3D

Keywords:
Tetrahedral mesh, Delaunay, Automatic mesher

Functional Description:
GHS3D is an automatic volume mesher
 URL:

Contact:
Paul Louis George

Participants:
Paul Louis George, Adrien Loseille, Frederic Alauzet
5.1.3 HEXOTIC

Keywords:
3D, Mesh generation, Meshing, Unstructured meshes, Octree/Quadtree, Multithreading, GPGPU, GPU

Functional Description:
Input: a triangulated surface mesh and an optional size map to control the size of inner elements.
Output: a fully hexahedral mesh (no hybrid elements), valid (no negative jacobian) and conformal (no dangling nodes) whose surface matches the input geometry.
The software is a simple command line that requires no knowledge on meshing. Its arguments are an input mesh and some optional parameters to control elements sizing, curvature and subdomains as well as some features like boundary layers generation.
 URL:

Contact:
Loic Marechal

Participant:
Loic Marechal

Partner:
Distene
5.1.4 Metrix

Name:
Metrix: Error Estimates and Mesh Control for Anisotropic Mesh Adaptation

Keywords:
Meshing, Metric, Metric fields

Functional Description:
Metrix is a software that provides by various ways metric to govern the mesh generation. Generally, these metrics are constructed from error estimates (a priori or a posteriori) applied to the numerical solution. Metrix computes metric fields from scalar solutions by means of several error estimates: interpolation error, isolines error estimate, interface error estimate and goal oriented error estimate. It also contains several modules that handle meshes and metrics. For instance, it extracts the metric associated with a given mesh and it performs some metric operations such as: metric gradation and metric intersection.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.5 VIZIR

Name:
Interactive visualization of hybrid, curbed and highorder mesh and solution

Keyword:
Mesh

Functional Description:
Vizir is a light, simple and interactive mesh visualization software, including: (i) A curved meshes visualizator: it handles high order elements and solutions, (ii) Hybrid elements mesh visualization (pyramids, prisms, hexahedra), (iii) Solutions visualization : clip planes, capping, isolines, isosurfaces.
 URL:
 Publication:

Contact:
Adrien Loseille

Participants:
Adrien Loseille, Rémi Feuillet, Matthieu Maunoury
5.1.6 Wolf

Keyword:
Scientific calculation

Functional Description:
Numerical solver for the Euler and compressible NavierStokes equations with turbulence modelling. ALE formulation for moving domains. Modules of interpolation, mesh optimisation and moving meshes. Wolf is written in C++, and may be later released as an opensource library. FELiScE was registered in July 2014 at the Agence pour la Protection des Programmes under the Inter Deposit Digital Number IDDN.FR.001.340034.000.S.P.2014.000.10000.
 URL:

Contact:
Frederic Alauzet

Participants:
Frederic Alauzet, Adrien Loseille, Rémi Feuillet, Lucille Marie Tenkes, Francesco Clerici
5.1.7 WolfBloom

Keyword:
Scientific calculation

Functional Description:
WolfBloom is a structured boundary layer mesh generator using a pushing approach. It start from an existing volume mesh and insert a structured boundary layer by pushing the volume mesh. The volume mesh deformation is solved with an elasticity analogy. Meshconnectivity optimizations are performed to control volume mesh element quality.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, David Marcum, Frederic Alauzet
5.1.8 WolfElast

Keyword:
Scientific calculation

Functional Description:
WolfElast is a linear elasticity solver using the P1 to P3 FiniteElement method. The Young and Poisson coefficient can be parametrized. The linear system is solved using the Conjugate Gradient method with the LUSGS preconditioner.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.9 WolfInterpol

Keyword:
Scientific calculation

Functional Description:
WolfInterpol is a tool to transfer scalar, vector and tensor fields from one mesh to another one. Polynomial interpolation (from order 2 to 4) or conservative interpolation operators can be used. WolfInterpol also extract solutions along lines or surfaces.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.10 WolfMovMsh

Keyword:
Scientific calculation

Functional Description:
WolfMovMsh is a moving mesh algorithm coupled with meshconnectivity optimization. Mesh deformation is computed by means of a linear elasticity solver or a RBF interpolation. Smoothing and swapping mesh optimization are performed to maintain good mesh quality. It handles rigid bodies or deformable bodies, and also rigid or deformable regions of the domain. Highorder meshes are also handled
 URL:

Contact:
Paul Louis George

Participants:
Adrien Loseille, Frederic Alauzet
5.1.11 WolfNsc

Keyword:
Scientific calculation

Functional Description:
WolfNsc is numerical flow solver solving steady or unsteady turbulent compressible Euler and NavierStokes equations. The available turbulent models are the SpalartAlmaras and the Menter SST komega. A mixed finite volume  finite element numerical method is used for the discretization. Second order spatial accuracy is reached thanks to MUSCL type methods. Explicit or implicit time integration are available. It also resolved dual (adjoint) problem and compute error estimate for mesh adaptation.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.12 WolfShrimp

Keyword:
Scientific calculation

Functional Description:
WolfShrimp is a generic mesh partitioner for parallel mesh generation and parallel computation. It can partition planar, surface (manifold and non manifold), and volume domain. Several partitioning methods are available: Hilbertbased, BFS, BFS with restart. It can work with or without weight function and can correct the partitions to have only one connected component.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.13 WolfSpyder

Keyword:
Scientific calculation

Functional Description:
WolfSpyder is a metricbased highorder mesh quality optimizer using vertex smoothing and edge/face swapping.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
5.1.14 WolfXfem

Keyword:
Scientific calculation

Functional Description:
WolfXfem is a tool providing the mesh of the intersection between a surface mesh and a volume mesh.
 URL:

Contact:
Frederic Alauzet

Participants:
Adrien Loseille, Frederic Alauzet
6 New results
6.1 Mesh glossary
Participants: Frédéric Alauzet [correspondant], Paul Louis George [correspondant], Adrien Loseille.
After the publication of the third volume on Meshing, Geometric Modeling and Numerical Simulation 3 was published in 2020 14, 31, 29, books also available in French 13, 30, 28, a glossary in French on mesh has been written 6. From A to Z, this glossary provides definitions and gives a number of comments. Moreover, subtilities are given probably not so easy to understand even by people familiar with the french dialect.
6.2 Numerical simulations on GPU with the GMlib v3.0 library
Participants: Loïc Maréchal [correspondant], Julien Vanharen.
The whole library was completely rewritten to implement an automatic finiteelement shader generation that converts a simple user source code into an OpenCL source that is in compiled on the GPU at run time. The library handles all meshing data structures, from file reading, renumbering and vectorizing for efficient access on the GPU, and transfer to the graphic card, all automatically and transparently. With this framework, the user can focus on the calculation part of the code, known as kernel, as all the rest is taken care of by the library. The OpenCL language was chosen as it is hardware agnostic and runs on any GPU (Intel, Nvidia and AMD) and can also use the multicore and vector capacities of modern CPUs.
Julien Vanharen developed a basic heat solver using the v3.0 as a test case so we could validate the software with various boundary conditions, calculation scheme, unstructured meshes and different memory access patterns with success. Even with basic calculation which does not stress the full GPU's power, we achieved two orders of magnitude greater speed against a single CPU core and one order of magnitude compared to a multithreaded implementation. As Julien moved to ONERA, we plan on setting up a collaboration between the two teams to implement more complex HPC codes.
6.3 High Order Meshing: from a straight mesh to a curved one
Participants: Loïc Maréchal [correspondant].
Works continued on P1toPk, a software that transform any first order hybrid mesh (triangles, quads, tets, pyramids, prisms and hexes) into a second order one while respecting a prescribe surface curvature. Efforts were made on boundary layers curving, which was challenging because jacobian validity is harder to guarantee as the elements get highly stretched, and a lot effort were also made to speed up the code by optimizing mathematical operations and parallelizing them. The code is now mature enough to be sent to industrial users for real life usage and we are waiting for valuable feedback in the present year.
6.4 Pixelexact rendering for highorder meshes and solutions
Participants: Adrien Loseille [correspondant], Rémi Feuillet, Matthieu Maunoury.
Classic visualization software like ParaView 40, TecPlot 41, FieldView 46, Ensight 39, Medit 26, Vizir (OpenGL legacy based version) 52, Gmsh 36, ... historically rely on the display of linear triangles with linear solutions on it. More precisely, each element of the mesh is divided into a set of elementary triangles. At each vertex of the elementary triangle is attached a value and an associated color. The value and the color inside the triangle is then deduced by a linear interpolation inside the triangle. With the increase of highorder methods and highorder meshes, these softwares adapted their technology by using subdivision methods. If a mesh has highorder elements, these elements are subdivided into a set of linear triangles in order to approximate the shape of the highorder element 82. Likewise, if a mesh has a highorder solution on it, each element is subdivided into smaller linear triangles in order to approximate the rendering of the highorder solution on it. The subdivision process can be really expensive if it is done in a naive way. For this reason, mesh adaptation procedures 69, 57, 58 are used to efficiently render highorder solutions and highorder elements using the standard linear rendering approaches. Even when optimized these approaches do have a huge RAM memory footprint as the subdivision is done on CPU in a preprocessing step. Also the adaptive subdivision process can be dependent on the palette (i.e. the range of values where the solution is studied) as the color only vary when the associated value is in this range. In this case, a change of palette inevitably imposes a new adaptation process. Finally, the use of a non conforming mesh adaptation can lead to a discontinuous rendering for a continuous solution.
Other approaches are specifically devoted to highorder solutions and are based on ray casting 63, 64, 66. The idea is for a given pixel, to find exactly its color. To do so, for each pixel, rays are cast from the position of the screen in the physical space and their intersection with the scene determines the color for the pixel. If highorder features are taken into account, it determines the color exactly for this pixel. However, this method is based on two nonlinear problems: the rootfinding problem and the inversion of the geometrical mapping. These problems are really costly and do not compete with the interactivity of the standard linear rendering methods even when these are called with a subdivision process unless they are done conjointly on the GPU. However, synchronization between GPU and OpenGL buffer are nontrivial combination.
The proposed method intends to be a good compromise between both methods. It does guarantee pixelexact rendering on linear elements without extra subdivision or ray casting and it keeps the interactivity of a classical method. Moreover, the subdivision of the curved entities is done on the fly on GPU which leaves the RAM memory footprint at the size of the loaded mesh.
We are developing a software, ViZiR 4, with exact non linear solution rendering to address the highorder visualization challenge 1. ViZiR 4 is bundled as a light, simple and interactive highorder meshes and solutions visualization software. It is based on OpenGL 4 core graphic pipeline. The use of OpenGL Shading Language (GLSL) allows to perform pixel exact rendering of high order solutions on straight elements and almost pixel exact rendering on curved elements (highorder meshes). ViZiR 4 enables the representation of high order meshes (up to degree 4) and high order solutions (up to degree 10) with pixel exact rendering. Furthermore, in comparison with standard rendering techniques based on legacy OpenGL, the use of OpenGL 4 core version improves the speed of rendering, reduces the memory footprint and increases the flexibility. Many postprocessing tools, such as picking, hidding surfaces, isolines, clipping, capping, are integrated to enable on the fly the analysis of the numerical results.
We added in ViZiR 4 the visualization of polygonal and polyhedral meshes. Such meshes offer flexibility as the number of vertices and faces are arbitrary. Dual meshes are examples of such meshes. Many visualization softwares do not handle polygons and when it is the case, the interactivity is often limited. We showed that our graphic pipeline can be used to render polygons. New keywords have been introduced in the libMeshb library to define these polygons and polyhedra. Furthermore, some functions have been introduced to ease the access of these data. New tessellation algorithms have been developed and do not add extra vertices in the tessellation as the aim is to minimize the number of triangles to maximize the rendering performances.
6.5 Highorder mesh generation
Participants: Frédéric Alauzet [correspondant], Adrien Loseille, Rémi Feuillet, Dave Marcum, Lucien Rochery.
For years, the resolution of numerical methods has consisted in solving Partial Derivative Equations by means of a piecewise linear representation of the physical phenomenon on linear meshes. This choice was merely driven by computational limitations. With the increase of the computational capabilities, it became possible to increase the polynomial order of the solution while keeping the mesh linear. This was motivated by the fact that even if the increase of the polynomial order requires more computational resources per iteration of the solver, it yields a faster convergence of the approximation error 3 81 and it enables to keep track of unsteady features for a longer time and with a coarser mesh than with a linear approximation of the solution. However, in 17, 45, it was theoretically shown that for elliptic problems the optimal convergence rate for a highorder method was obtained with a curved boundary of the same order and in 12, evidence was given that without a highorder representation of the boundary the studied physical phenomenon was not exactly solved using a highorder method. In 85, it was even highlighted that, in some cases, the order of the mesh should be of a higher degree than the one of the solver. In other words, if the used mesh is not a highorder mesh, then the obtained highorder solution will never reliably represent the physical phenomenon.
Based on these issues, the development of highorder mesh generation procedures appears mandatory 22. To generate highorder meshes, several approaches exist. The first approach was tackled twenty years ago 20 for both surface and volume meshing. At this moment the idea was to use all the meshing tools to get a valid highorder mesh. The same problem was revisited a few years later in 75 for biomedical applications. In these first approaches and in all the following, the underlying idea is to use a linear mesh and elevate it to the desired order. Some make use of a PDE or variational approach to do so 7, 67, 23, 60, 80, 83, 37, others are based on optimization and smoothing operations and start from a linear mesh with a constrained highorder curved boundary in order to generate a suitable highorder mesh 44, 27, 78. Also, when dealing with NavierStokes equations, the question of generating curved boundary layer meshes (also called viscous meshes) appears. Most of the time, dedicated approaches are setup to deal with this problem 61, 43. In all these techniques, the key feature is to find the best deformation to be applied to the linear mesh and to optimize it. The prerequisite of these methods is that the initial boundary is curved and will be used as an input data. A natural question is consequently to study an optimal position of the highorder nodes on the curved boundary starting from an initial linear or highorder boundary mesh. This can be done in a coupled way with the volume 70, 79 or in a preprocessing phase 71, 72. In this process, the position of the nodes is set by projection onto the CAD geometry or by minimization of an error between the surface mesh and the CAD surface. Note that the vertices of the boundary mesh can move as well during the process. In the case of an initial linear boundary mesh with absence of a CAD geometry, some approaches based on normal reconstructions can be used to create a surrogate for the CAD model 82, 38. Finally, a last question remains when dealing with such highorder meshes: Given a set of degrees of freedom, is the definition of these objects always valid ? Until the work presented in 33, 42, 32, no real approach was proposed to deal in a robust way with the validity of highorder elements. The novelty of these approaches was to see the geometrical elements and their Jacobian as Bézier entities. Based on the properties of the Bézier representation, the validity of the element is concluded in a robust sense, while the other methods were only using a sampling of the Jacobian to conclude about its sign without any warranty on the whole validity of the elements.
In this context, several issues have been addressed : the analogy between highorder and Bézier elements, the development of highorder error estimates suitable for parametric highorder surface mesh generation and the generalization of mesh optimization operators and their applications to curved mesh generation, movingmesh methods, boundary layer mesh generation and mesh adaptation.
Metric fields are the link between particular error estimates  be they for loworder 47, 48 or highorder methods 19, for the solution of a PDE 50 or a quantity of interest derived from it such as drag or lift 51  and automatic mesh adaptation. In the case of linear meshes, a metric field locally distorts the measure or distance such that, when the mesh adaptation algorithm has constructed an uniform mesh in the induced Riemannian space, it is strongly anisotropic in the usual Euclidean (physicial) space. As such, anisotropy arises naturally, without it ever being explicitely sought by the (re)meshing algorithm.
We seek to extend these principles of metricbased ${P}^{1}$ adaptation to highorder meshes. In particular, we expect the meshing process to naturally recover curvature from the variations of the metric field, very much like ${P}^{1}$ remeshing recovers anisotropy from local values of the metric field. As such, curvature must be the consequence of a simple geometric property computed in the Riemannian space, like anisotropy is the consequence of unitness in a space where distances are distorted. Therefore, we propose Riemannian edge length minimization (or geodesic seeking as in 84) as the driver for metric field curvature recovery.
The metric field's own intrinsic curvature may derive from any error estimate, be it boundary approximation error 24, 21 or an interpolation error estimate. So far, interpolation error estimates on highorder elements are limited to isotropy (16 in ${L}^{2}$ and 62 in ${L}^{1}$ norms) or require that the curvature of the element be bounded, essentially establishing a range where it may be considered linear 10.
Robustness and modularity of the general remeshing algorithm may be derived from the use of a single topological operator such as the cavity operator 49, 53, 54. This operator remeshes element subsets (socalled cavities) by reconnecting cavity boundary nodes to a given point in space, already belonging to the mesh or otherwise. This very elementary operation can handle the most common topological operations: insertions, collapses, edge or face swaps. Therefore, it is central both to mesh adaptation (node creation, deletion) and to mesh optimization (mainly swaps).
A new ${P}^{2}$ cavity operator was devised and implemented in $3D$. Volume curvature is driven, first, by a purely metricbased curving procedure using a Riemannian edge length minimization algorithm. The considered costfunction is strongly nonlinear, as it depends on the geometry of the edge as well as on the metric field along the edge which, itself, depends on the geometry of the edge. Analytical derivatives of this quantity with regards to edge shape were computed nonetheless, which was necessary to implement optimization in $3D$ with reasonable CPU costs. This first required the extension of the logeuclidean metric interpolation scheme to highorder elements. A version for simplices of all dimensions and degrees was proposed based on Lagrange interpolation of the logmetrics. This means that the metric field along the ${P}^{2}$ edge can now be described in terms of the metrics at the two extremities and at the Lagrange control point. The dependency of the metric at the Lagrange control point must in turn be identified. It is based on localization of the point in a socalled backmesh and, again, logeuclidean metric interpolation in the host back element. At last, derivatives of Riemannian edge length could be computed, and were used to implement a Newton minimization algorithm that is capable of finding local minima of edge length at a rate of 4e4 edges per second while remaining entirely consistent with the logeuclidean framework used in our anisotropic remeshing algorithm. Secondly, final cavities are made valid or improved using a new Jacobian correction algorithm that can very quickly optimize entire shells of edges using a recently identified property of Jacobian determinant control coefficients, namely that they depend linearly on a given control point. Using the very powerful simplex algorithm for linear programs, we can find a global solution to maximizing the minimum Jacobian of all elements sharing an edge in a single pass (with average runtime of 5e5 seconds on a laptop). Three communications have been made on the subject in the last year, with peerreviewed proceedings 3, 4, 2. Applications on large meshes deriving from industrial cases were carried out, illustrating the ability of these optimization procedures to curved meshes of 20M elements in under 10min.
6.6 Unstructured anisotropic mesh adaptation for 3D RANS turbomachinery applications
Participants: Frédéric Alauzet [correspondant], Adrien Loseille [correspondant], Julien Vanharen.
The scope of this paper is to demonstrate the viability and efficiency of unstructured anisotropic mesh adaptation techniques to turbomachinery applications. The main difficulty in turbomachinery is the periodicity of the domain that must be taken into account inthe solution meshadaptive process. The periodicity is strongly enforced in the flow solver using ghost cells to minimize the impact on the source code. For the mesh adaptation, the local remeshing is done in two steps. First, the inner domain is remeshed with frozen periodic frontiers, and, second, the periodic surfaces are remeshed after moving geometrice ntities from one side of the domain to the other. One of the main goal of this work is to demonstrate how mesh adaptation, thanks to its automation, is able to generate meshes that are extremely difficult to envision and almost impossible to generate manually. This study only considers featurebased error estimate based on the standard multiscale Lpinterpolation error estimate. We presents all the specific modifications that have been introduced in the adaptive process to deal with periodic simulations used for turbomachinery applications. The periodic mesh adaptation strategy is then tested and validated on the LS89 high pressure axial turbine vane and the NASA Rotor 37 test cases.
6.7 Mixedelement mesh adaptation for CFD simulations
Participants: Frédéric Alauzet [correspondant], Lucille Marie Tenkès, Julien Vanharen, Adrien Loseille, Cosimo Tarsia Morisco.
Due to their various nature, physical phenomena that we seek to capture in CFD simulations may have specific specific mesh requirements. For example, to solve the boundary layer, some numerical schemes favor structured meshes respecting alignment with the boundary of the domain, while these constraints are not necessary elsewhere. Our approach is to use the techniques of metricbased mesh adaptation to generate a mixedelement mesh that can fulfill these different mesh requirements. This approach is based on the metricorthogonal pointplacement, creating some structured parts from the intrinsic directional information bore by the metricfield. Some unstructured areas may remain where structure is not needed. The main goals of this work are to improve the orthogonality of the output mesh and its alignment with the metric field. This work has three main axes. First, we have improved the preprocessing gradation step to smooth the metric field and improve the orthogonality of the final mesh. Then, we have studied two methods to obtain quadrilaterals: one using an a priori quadrilaterals recombination, the other detecting straightforwardly the orthogonal patterns during the remeshing step. Finally, the work on the solver Wolf has been carried on and corrected to perform robust and accurate simulations on mixedelement meshes. These three developments were embodied in a mixedelement adaptation loop. The first two topics are detailed in what follows.
Enhanced metric gradation correction
The previously described generation method highly relies on the metric field. However, a metric field computed from a solution during the adaptation process is most of the time quite messy and shows abrupt size variations. In standard mesh adaptation, it leads to lowquality elements. In orthogonal mesh adaptation, it additionally breaks the alignment and the structure of the output mesh. An additional step to smooth the input metric field is therefore required. In the context of mixedelement mesh adaptation, this gradation correction process has been modified to improve the number and the quality of the quadrilaterals in the final mesh. Further developments have been considered on this topic, in particular to increase the robustness of the method. Results have been published in 77.
A posteriori and a priori mesh generation
Metricorthogonal pointplacement is currently used to generate quasistructured meshes with rightangled triangles where the metric is the most anisotropic and unit triangles elsewhere. The aim of this work is to recover some quadrilaterals in the structure. To do so, two approaches can be considered: an a posteriori quadrilateral recombination based on geometrical criteria, and an a priori quadrilateral detection. The latter is more straightforward because it uses directly the pointplacement information. A framework was established to set up this method. Developments and preliminary results were presented in 5.
In order to obtain a correct metric field on hybrid meshes, a robust hybrid solver is mandatory. When dealing with 2D (3D) elements different from triangles (tetrahedra), the most tricky aspect is the gradient formulation. This is due to the fact that within a Finite Elements interpolation framework, the gradient on an element with more than three nodes (i.e. not simplicial complex) is not elementwise constant. This brings many added difficulties to the flux balance computation. A first attempt at performing inviscid and laminar simulations on hybrid meshes was to approximate gradients on quadrilaterals with its isobarycenter values 5. The extension of this formulation to turbulent flows, however highlighted a lack of efficiency and robustness. For this reason, a APFE (APproximated FiniteElement) method 68 has been implemented and as well extended to a implicit time integration scheme. This approach turns out to be very efficient and robust in many fullystructured mesh verification cases. The extension to 3D cases (prisms and pyramids) is ongoing.
6.8 Coupled flow and adjoint solver
Participants: Francesco Clerici [correspondant], Frédéric Alauzet.
When solving the RANS equations, usually one decouples the equations relative to the meanflow and the equations relative to turbulence. This division provides two separated systems to be solved at each time step, one relative to the meanflow and the other relative to the turbulence. This presents two main drawbacks: in the flow solver, the Jacobian of the system lacks of the terms bounding the meanflow and the turbulence, and this can slow down the residual convergence. The second drawback regards the adjoint problem, which consists into a linear system assembled with the transpose of the Jacobian matrix of the residuals and, on the righthand side, the derivative of an aeronautical coefficient with respect to the flow variables. A Jacobian missing the coupling terms between the meanflow and the turbulence provides a null adjoint turbulent viscosity, and this is a limitation in the development of more complex discretization error estimates. We have therefore developed a 2D version of the coupled flow and adjoint solver which includes in the Jacobian the coupling terms between the meanflow and the turbulence. When have tested this method on the 2D geometry provided for the 4th CFD AIAA High Lift Prediction Workshop, and the result provided several features to emerge, such as a high mesh refinement inside the boundary layer of the leading edges, and inside the regions of high turbulence destruction. The work is pursued in collaboration with Philippe Spalart (Boeing), and is presented in 18.
6.9 Turbulent error estimate
Participants: Francesco Clerici [correspondant], Frédéric Alauzet.
Goaloriented mesh adaptation is a methodology used to adapt the mesh in order to minimize the discretization error commited on a functional depending on the solution. As an intermediate step, one finds an upper bound to such discretization error taking the form of a weighted sum of the interpolation errors of the solution, and such upper bound is called error estimate. Regarding Wolf and the RANS equations, up to now we have focused only on the meanflow part of such an error estimate, meaning that the terms coming from turbulence have been neglected. The scope of this work is to enrich the error estimate with the information coming from the turbulence. In particular, the methodology has been tested on the 2D geometry provided for the 4th CFD AIAA High Lift Prediction Workshop, providing high mesh refinements on the boundaries of the turbulent regions. Also this work is pursued in collaboration with Philippe Spalart (Boeing), and is presented in 18.
7 Bilateral contracts and grants with industry
7.1 Bilateral contracts with industry
Participants: Frédéric Alauzet [correspondant], Adrien Loseille.
 Boeing
 Safran Tech
7.2 Bilateral grants with industry
Participants: Frédéric Alauzet [correspondant], Adrien Loseille.
 Projet RAPID DGA
ANR IMPACTS 20182021
Ideal Mesh generation for modern solvers and comPuting ArchiteCTureS.
 Coordinateur : Adrien Loseille
 The rapid improvement of computer hardware and physical simulation capabilities has revolutionized science and engineering, placing computational simulation on an equal footing with theoretical analysis and physical experimentation. This rapidly increasing reliance on the predictive capabilities has created the need for rigorous control of numerical errors which strongly impact these predictions. A rigorous control of the numerical error can be only achieved through mesh adaptivity. In this context, the role of mesh adaptation is prominent, as the quality of the mesh, its refinement, and its alignment with the physics are major contributions to these numerical errors. The IMPACTS project aims at pushing the envelope in mesh adaptation in the context of large size, very high fidelity simulations by proposing a new adaptive mesh generation framework. This framework will be based on new theoretical developments on Riemannian metricfield and on innovative algorithmic developments coupling a unique cavityoperator with an advancingpoint techniques in order to produce high quality hybrid, curved and adapted meshes.
8 Scientific production
8.1 Publications of the year
International journals
 1 articleOn pixelexact rendering for highorder mesh and solution.Journal of Computational Physics424January 2021, 109860
International peerreviewed conferences
 2 inproceedingsDevelopments on the P2 cavity operator and Bézier Jacobian correction using the simplex algorithm..AIAA SCITECH 2022 ForumSan Diego, United StatesAmerican Institute of Aeronautics and AstronauticsJanuary 2021
 3 inproceedingsP2 Cavity Operator and Riemannian Curved Edge Length Optimization: a Path to HighOrder Mesh Adaptation.AIAA Scitech Forum 2021Virtual, United StatesAmerican Institute of Aeronautics and AstronauticsJanuary 2021
 4 inproceedingsP2 Cavity Operator with MetricBased Volume and Surface Curvature.29th International Meshing Roundtable (IMR)Virtual, FranceJune 2021
 5 inproceedingsHybrid anisotropic mesh adaptation using metricorthogonal approach.AIAA Scitech 2021 ForumVirtual, United StatesAmerican Institute of Aeronautics and AstronauticsJanuary 2021
Reports & preprints
 6 reportLe glossaire du maillage.RT0515INRIASeptember 2021
8.2 Cited publications
 7 articleA method for computing curved meshes via the linear elasticity analogy, application to fluid dynamics problems.International Journal for Numerical Methods in Fluids7642014, 246266

8
articleA decade of progress on anisotropic mesh adaptation for Computational Fluid Dynamics.
722016, 1339 
9
articleHigh Order Sonic Boom Modeling by Adaptive Methods.
2292010, 561593  10 bookAnisotropic finite elements: local estimates and applications.3Teubner Stuttgart1999

11
articleMetricbased anisotropic mesh adaptation for threedimensional timedependent problems involving moving geometries.
3312017, 157187  12 articleHighorder accurate discontinuous finite element solution of the 2D Euler equations.Journal of computational physics13821997, 251285
 13 bookMaillage, modélisation géométrique et simulation numérique 1: Fonctions de forme, triangulations et modélisation géométrique.1ISTE Group2017
 14 bookMeshing, Geometric Modeling and Numerical Simulation 1: Form Functions, Triangulations and Geometric Modeling.John Wiley & Sons2017
 15 articleParametric surface meshing using a combined advancingfront  generalizedDelaunay approach.International Journal for Numerical Methods in Engineering49122000, 233259
 16 articleInfluence of ReferencetoPhysical Frame Mappings on Approximation Properties of Discontinuous Piecewise Polynomial Spaces.Journal of Scientific Computing5209 2012
 17 incollectionThe combined effect of curved boundaries and numerical integration in isoparametric finite element methods.The mathematical foundations of the finite element method with applications to partial differential equationsElsevier1972, 409474
 18 articleCoupled adjoint solver and turbulent error estimate for anisotropic mesh adaptation in highfidelity RANS simulations..AIAA J.2022
 19 articleVery High Order Anisotropic MetricBased Mesh Adaptation in 3D.Procedia Engineering16325th International Meshing Roundtable2016, 353  365URL: http://www.sciencedirect.com/science/article/pii/S187770581633380X
 20 inproceedingsCurvilinear Mesh Generation in 3D.Proceedings of the 7th International Meshing Roundtable1999, 407417
 21 inproceedingsAnisotropic Error Estimate for Highorder Parametric Surface Mesh Generation.28th International Meshing RoundtableBuffalo, NY, United StatesOctober 2019
 22 articleOptimization of P2 meshes and applications.ComputerAided Design124April 2020, 102846
 23 articleHighorder Unstructured Curved Mesh Generation Using the Winslow Equations.J. Comput. Phys.307February 2016, 114
 24 miscAbout Surface Remeshing.2000
 25 inproceedingsAbout surface remeshing.Proceedings of the 9th international meshing roundtableNew Orleans, LO, USA2000, 123136
 26 miscMedit: An interactive mesh visualization software, INRIA Technical Report RT0253.2001
 27 incollectionDefining quality measures for mesh optimization on parameterized CAD surfaces.Proceedings of the 21st International Meshing RoundtableSpringer2013, 85102
 28 bookMaillage, modélisation géométrique et simulation numérique 3: Stockage, transformation, utilisation et visualisation de maillage.4ISTE Group2020
 29 bookMeshing, Geometric Modeling and Numerical Simulation 3: Storage, Visualization and In Memory Strategies.John Wiley & Sons2020
 30 bookMaillage, modélisation géométrique et simulation numérique 2: Métriques, maillages et adaptation de maillages.2ISTE Group2018
 31 bookMeshing, Geometric Modeling and Numerical Simulation, Volume 2: Metrics, Meshes and Mesh Adaptation.John Wiley & Sons2019
 32 articleGeometric validity (positive jacobian) of highorder Lagrange finite elements, theory and practical guidance.Engineering with computers3232016, 405424
 33 articleConstruction of tetrahedral meshes of degree two.International Journal for Numerical Methods in Engineering9092012, 1156,1182

34
article``Ultimate'' robustness in meshing an arbitrary polyhedron.
5872003, 10611089 
35
articleAutomatic mesh generator with specified boundary.
921991, 269288  36 articleGmsh: A 3D finite element mesh generator with builtin pre and postprocessing facilities.International Journal for Numerical Methods in Engineering79112009, 13091331
 37 articleGeneration of unstructured curvilinear grids and highorder discontinuous Galerkin discretization applied to a 3D highlift configuration.International Journal for Numerical Methods in Fluids8262016, 316333
 38 articleAutomated loworder to highorder mesh conversion.Engineering with Computers351Jan 2019, 323335
 39 miscEnsight.
 40 miscParaView.
 41 miscTecPlot.
 42 articleGeometrical validity of curvilinear finite elements.Journal of Computational Physics233152013, 359372
 43 inproceedingsCurving for Viscous Meshes.27th International Meshing RoundtableChamSpringer International Publishing2019, 303325
 44 incollectionHighOrder Mesh Curving Using WCN Mesh Optimization.46th AIAA Fluid Dynamics ConferenceAIAA AVIATION ForumAmerican Institute of Aeronautics and Astronautics2016
 45 articleOptimal isoparametric finite elements and error estimates for domains involving curved boundaries.SIAM journal on numerical analysis2331986, 562580
 46 miscFieldView.
 47 articleContinuous mesh framework part I: wellposed continuous interpolation error.SIAM Journal on Numerical Analysis4912011, 3860
 48 articleContinuous mesh framework part II: validations and applications.SIAM Journal on Numerical Analysis4912011, 6186
 49 articleUnique cavitybased operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes.ComputerAided Design852017, 5367
 50 phdthesisAnisotropic 3D hessianbased multiscale and adjointbased mesh adaptation for Computational fluid dynamics: Application to high fidelity sonic boom prediction..Université Pierre et Marie Curie  Paris VIDecember 2008
 51 articleFully anisotropic goaloriented mesh adaptation for 3D steady Euler equations.Journal of Computational Physics22982010, 2866  2897URL: http://www.sciencedirect.com/science/article/pii/S0021999109007001
 52 miscAn introduction to Vizir: an interactive mesh visualization and modification software.2016
 53 inbookCavityBased Operators for Mesh Adaptation.51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace ExpositionURL: https://arc.aiaa.org/doi/abs/10.2514/6.2013152
 54 inbookRecent Improvements on CavityBased Operators for RANS Mesh Adaptation.2018 AIAA Aerospace Sciences MeetingURL: https://arc.aiaa.org/doi/abs/10.2514/6.20180922
 55 inproceedingsA new approach to octreebased hexahedral meshing.2001, 209221

56
inproceedingsAdvances in OctreeBased AllHexahedral Mesh Generation: Handling Sharp Features.
18Salt Lake City, UT, USA2009, 6584  57 articleWellsuited and adaptive postprocessing for the visualization of hp simulation results.Journal of Computational Physics3752018, 1179  1204
 58 phdthesisMéthode de visualisation adaptée aux simulations d'ordre élevé. Application à la compressionreconstruction de champs rayonnés pour des ondes harmoniques..Université de ToulouseFebruary 2019

59
inproceedingsResults from the 3rd Drag Prediction Workshop using NSU3D unstructured mesh solver.
45AIAA20070256, Reno, NV, USAJan 2007  60 articleHighorder curvilinear meshing using a thermoelastic analogy.ComputerAided Design722016, 130  139
 61 articleAn isoparametric approach to highorder curvilinear boundarylayer meshing.Computer Methods in Applied Mechanics and Engineering2832015, 636  650
 62 articleInterpolation error bounds for curvilinear finite elements and their implications on adaptive mesh refinement.Journal of Scientific Computing7822019, 10451062
 63 articleGPUBased Interactive CutSurface Extraction From HighOrder Finite Element Fields.IEEE Transactions on Visualization and Computer Graphics17122011, 180311
 64 articleElVis: A System for the Accurate and Interactive Visualization of HighOrder Finite Element Solutions.IEEE Transactions on Visualization and Computer Graphics18122012, 23252334

65
inproceedingsUnstructured Grid Adaptation: Status, Potential Impacts, and Recommended Investments Toward CFD Vision 2030.
46 20163323, Washington, D.C., USA2016  66 incollectionHighOrder Visualization with ElVis.Notes on Numerical Fluid Mechanics and Multidisciplinary DesignSpringer International Publishing2015, 521534
 67 inproceedingsCurved mesh generation and mesh refinement using Lagrangian solid mechanics.47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition2009, 949
 68 articleDiscretisation of diffusive fluxes on hybrid grids.Journal of Computational Physics22952010, 14251447
 69 articleEfficient visualization of highorder finite elements.International Journal for Numerical Methods in Engineering6952007, 750771
 70 articleHighorder mesh curving by distortion minimization with boundary nodes free to slide on a 3D CAD representation.ComputerAided Design7223rd International Meshing Roundtable Special Issue: Advances in Mesh Generation2016, 52  64

71
articleDefining an
${L}^{2}$ disparity Measure to Check and Improve the Geometric Accuracy of Noninterpolating Curved Highorder Meshes.Procedia Engineering1242015, 122134  72 articleGeneration of curved highorder meshes with optimal quality and geometric accuracy.Procedia engineering1632016, 315327
 73 articleSummary of the first AIAA CFD HighLift Prediction Workshop.Journal of Aircraft4862011, 20682079
 74 articleOverview and Summary of the Second AIAA High Lift Prediction Workshop.Journal of Aircraft5242015, 10061025
 75 articleMesh generation in curvilinear domains using highorder elements.International Journal for Numerical Methods in Engineering5312002, 207223
 76 techreportCFD Vision 2030 Study: A path to revolutionary computational aerosciences.NASAMarch 2014
 77 incollectionSize Gradation Control for Anisotropic Hybrid Meshes.Numerical Geometry, Grid Generation and Scientific Computing143Lecture Notes in Computational Science and EngineeringSpringer International PublishingMay 2021, 127139
 78 articleRobust untangling of curvilinear meshes.Journal of Computational Physics2542013, 8  26
 79 articleOptimizing the geometrical accuracy of curvilinear meshes.Journal of Computational Physics3102016, 361  380
 80 articleA Variational Framework for Highorder Mesh Generation.Procedia Engineering163Supplement C25th International Meshing Roundtable2016, 340  352
 81 phdthesisHighorder numerical methods for unsteady flows around complex geometries.Université de Toulouse2017
 82 inproceedingsCurved PN Triangles.Proceedings of the 2001 Symposium on Interactive 3D Graphics2001, 159166
 83 articleThe generation of arbitrary order curved meshes for 3D finite element analysis.Computational Mechanics513Mar 2013, 361374
 84 inproceedingsCurvilinear mesh adaptation.International Meshing RoundtableSpringer2018, 5769
 85 inproceedingsOn the Necessity of Superparametric Geometry Representation for Discontinuous Galerkin Methods on Domains with Curved Boundaries.23rd AIAA Computational Fluid Dynamics Conference2017