text:
Stéphanie Péron
email: stephanie.peron@onera.fr

Dylan Jude

Andrew Wissink

Ronan Boisard

St Éphanie P Éron

Assessment of reduced order fuselage and blade models for rotorcraft international aerodynamics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1.INTRODUCTION

Prediction of aircraft aerodynamics is a key component in the rotorcraft design process. Rotorcraft comprehensive analysis tools, such as RCAS [START_REF] Saberi | RCAS Overview of Capabilities, Validations, and Applications to Rotorcraft Problems[END_REF] and HOST, [START_REF] Benoit | HOST, a general helicopter simulation tool for Germany and France[END_REF] are tightly integrated into the design process of a new rotorcraft. However their aerodynamic models for the rotor and fuselage are often lowfidelity. For traditional, single-main-rotor helicopters, a low-fidelity model tuned for certain flight conditions may produce adequate results for the purposes of design. Modern rotorcraft configurations, however, often include multiple rotors and lifting surfaces with non-trivial interactional aerodynamics and low-fidelity models often cannot account for these factors. For the design process to include effects of complex interactional aero-dynamics, higher-fidelity aerodynamic models need to be introduced into the design loop.

At the opposite end of the spectrum from low-fidelity aerodynamic lies high-fidelity, computationally expensive tools using overset Computational Fluid Dynamics (CFD). Codes such as elsA [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF] (ONERA) and Helios [START_REF] Wissink | New Capabilities in CREATE-AV Helios Version 11[END_REF][START_REF] Jude | An octree-based, Cartesian Navier-Stokes solver for modern cluster architectures[END_REF] (CREATE TM -AV) are established, accurate CFD frameworks for simulating complex aerodynamics. Both frameworks use a "Chimera" approach by which unstructured, structured curvilinear, or strand near-body meshes are overset within a Cartesian background mesh. The near-body meshes resolve the near-wall viscous boundary-layer flow and the off-body Cartesian region prescribes Adaptive Mesh Refinement (AMR) to resolve aircraft or rotor wakes and inter-actional aerodynamics. A Helios or elsA simulation of an aircraft will commonly have tens or hundreds-of-millions of degrees of freedom, requiring parallel execution on a large scale high performance computing (HPC) supercomputer. Even with thousands of CPU cores on modern supercomputers, complex simulations can still take days or weeks to complete.

Research groups at ONERA and the US Army have recently developed mid-fidelity tools intended to fill the gap between the traditional low and high-fidelity tools available today. Specifically, these tools are designed to give reasonably accurate aerodynamic loads analysis of the interactional aerodynamics experienced by modern rotorcraft vehicles but at a much faster turnaround rate than high-fidelity overset body-fitted CFD methods. In both models, complex geometries that would normally require large near-body grids are replaced with an immersed boundary representation. The Immersed Boundary Method (IBM) interacts directly with the background Cartesian mesh to inject influence from body so that interactional effects are captured. Since viscous body-fitted meshes are not required, solutions are performed only on computationallyefficient Cartesian grids. Further, larger timesteps are allowed through elimination of cell size limitations of viscous boundary layer meshes. The combination of these two features enable accelerations of up to 40X-50X faster over traditional high-fidelity overset CFD methods.

The IBM approach does not work well for modeling rotor blades or other lifting surfaces without very fine Cartesian refinement. Helios solves this problem by representing each blade as an actuator line. Sectional loads based on 2D aerodynamics are integrated along the blade span and injected into the background Cartesian mesh as source terms in the governing Navier-Stokes equations.

The IBM and actuator line methods in ONERA considers a different approach within FAST/Cassiopee framework: URANS solution is computed on wall-resolved curvilinear grids describing the blades, in relative motion with respect to the off-body Cartesian mesh, in which the fuselage is modeled by IBM. A wall function is applied to recover the skin friction onto the fuselage, as detailed in P éron et al. [START_REF] P Éron | An immersed boundary method on Cartesian adaptive grids for the simulation of compressible flows around arbitrary geometries[END_REF] and Constant et al. [START_REF] Constant | An improved Immersed Boundary Method for turbulent flow simulations on Cartesian grids[END_REF] . The present study investigates the usage of Helios/ROAM and FAST/IBM to assess interactional aerodynamics of the 7A four-bladed rotor on the Dauphin fuselage (D365N test case) for which existing body-fitted CFD results are available 8 .

2.DESCRIPTION OF THE TEST-CASE

A scaled model of the Dauphin helicopter equipped with a powered main rotor was designed and tested in the late 1980s in the ONERA S2Ch wind tunnel and in ONERA F1 subsonic wind tunnel. The test provided experimental data, including steady and unsteady fuselage surface pressures, balance measurements, 2D and 3D PIV wake measurements. The tested helicopter model is a 1/7.7 Dauphin model with a 1.5 meter fuselage. The powered main rotor is a 0.75 meter radius 4-bladed rotor articulated in pitch, flap and lead-lag motions. The 0.05 m constant chord blade uses OA209 airfoils and has -12 deg/R of linear twist. The blade root cutout is at 0.275R. The clockwise rotor rotation is ensured by an electric engine which allows a blade tip speed of 100 m.s ´1 (1272 RPM). The F1 experimental fuselage model varies from the theoretical Dauphin 365N model definition, especially in the area of the engine inlet and exhaust fairings. Neither the tail rotor nor the rotor hub are included in the present simulations; the strut and wind tunnel walls have also been excluded.

In the F1 experiments, the D365N configuration is trimmed to a specific thrust (Z b =14.5), axial force (X b =0) and rotor lateral flapping β 1s =0. The baseline flow conditions are V 8 =15 m.s ´1, corresponding to a Mach number M a 8 =0.044 and a moderate advance ratio µ " 0.15, ´3˝f uselage angle of attack and zero sideslip, well suited for the study of the rotor-fuselage interaction, as the rotor wake impinges the fuselage and the empennage. The Reynolds number is 1.07 million per meter. The rotor control angles (collective, flap, lag and harmonics) are obtained from the rotorcraft comprehensive code HOST [START_REF] Benoit | HOST, a general helicopter simulation tool for Germany and France[END_REF] . Rigid blade deflections are applied about lag and flap hinges using three harmonics.

3.COMPUTATIONAL FRAMEWORK

3.1.Helios ROAM

Helios 4 is a high-fidelity coupled CFD/CSD infrastructure dedicated to aeromechanics predictions of rotorcraft. It relies traditionally on body conforming unstructured or structured-curvilinear meshes in the nearbody region and Cartesian meshes in the off-body region. The near-body grids rotate, move, and deform with the rotor blades and resolve the unsteady turbulent viscous scales in the boundary layer. The stationary off-body Cartesian solver resolves the wake through a combination of AMR and highorder accurate discretizations. Data transfer between the two mesh systems are performed through the use of overset domain connectivity transfer technology that dynamically cuts holes and fringes in the meshes to accommodate the relative motion between the mesh systems. Typical simulations may require several days on HPC systems while larger simulations requiring many rotor revolutions may take over a week to complete.

A medium-fidelity option has recently been implemented in Helios, named ROAM (Reduced-Order Aerodynamics Model), relying on Cartesian grids [START_REF] Wissink | New Capabilities in CREATE-AV Helios Version 11[END_REF] . It applies an actuator line model for the rotor and an immersed boundary representation of the fuselage. The Cartesian solver used by ROAM is the same off-body solver used to resolve the wake in near-body/off-body overset calculations in Helios. Calculations with ROAM are significantly faster than the traditional body-fitted calculations for two reasons; 1) it avoids the cost of near-body viscous boundary layer solution and domain connectivity because it does not explicitly resolve this region; and 2) it allows for much larger timesteps in the unsteady simulation. The methodology for using actuator lines in ROAM is based on lifting-line theory and similar to the approaches from by Mittal et al. [START_REF] Mittal | Improvements to the Actuator Line Modeling for Wind Turbines[END_REF] and Churchfield et al. [START_REF] Churchfield | An Advanced Actuator Line Method for Wind Energy Applications and Beyond[END_REF] . A lifting line along the blade is discretized into a set of control points and the flow induced velocities at these control points are interpolated from the Cartesian mesh system modeling the flow field. To exclude the effects of bound circulation, a circle of sensors is created around each control point and the velocity at the control point is found as an average of all sensors. Once the velocities at control points are obtained, linear aerodynamic theory augmented with airfoil tables are used to determine the force coefficients and compute the force distribution. The force distribution is then embedded into the flow field equations as source terms for the momentum and energy equations. An illustration of the lifting line in a background mesh is shown in Fig. 1. The normal component of the resultant force is injected back into the domain at both the 1/4 and 3/4 chord locations in proportions to match the aerodynamic pitching moment at that span-wise section.

Helios has two background Cartesian AMR solvers available for use. SAMCart uses patch-based refinement and Orchard uses an octree-based data structure for refinement. The underlying numerics for the Cartesian Navier-Stokes equation are the same between both approaches, however the octree-based approach in Orchard results in performance and scalability improvements over SAMCart. For all Helios cases in this work Orchard is used as the background Cartesian solver. The inviscid terms of the Navier-Stokes equations use a 5 thorder central stencil with scalar dissipation and the viscous terms use 4 th -order stencils. The temporal order of accuracy is second order using an implicit backwards differencing algorithm with dual-timestepping. Some high-fidelity results in this work are obtained for reference from Helios using body-fitted and the mStrand [START_REF] Lakshminarayan | Development and validation of a multi-strand solver for complex aerodynamic flows[END_REF] solver overset in a similar background domain.

3.2.FAST/IBM

FAST solver 12 is a structured multi-block Finite-Volume solver devoted to high-fidelity simulations (including DNS, LES and RANS-LES simulations), developed by ONERA. Current FAST architecture relies on an hybrid MPI/OpenMP parallelism on CPUs and is integrated within a CGNS/Python framework in combination with Cassiopee set of modules for pre, co-and post-processing [START_REF] Benoit | Cassiopee: A CFD pre-and postprocessing tool[END_REF] , providing a versatile simulation tool to perform either high-fidelity unsteady simulations on structured multi-block meshes, Cartesian IBM meshes, and more recently hybrid IBM and body-fitted meshes together. FAST Cartesian solver is applied on Carte-sian grids, dividing by two the memory and CPU time compared to the curvilinear solver.

An overset grid capability for rotorcraft applications has been recently developed within Cassiopee and applied with FAST solver. For that range of applications, curvilinear grids are defined around rotor blades whereas the off-body region is meshed by an octree-like Cartesian mesh, as described in P éron et al. [START_REF] P Éron | Automatic off-body overset adaptive Cartesian mesh method based on an octree approach[END_REF] . For that purpose, Chimera transfers are performed at each sub-time step between curvilinear grids and Cartesian grids, hole-cutting and data transfers between overset grids have been especially designed for this type of mesh, where the Cartesian mesh is dominant and where curvilinear grids represent bodies in relative motion.

Here, unsteady Reynolds-Averaged Navier-Stokes equations (URANS) are solved, using Spalart-Allmaras turbulence model. 2 nd -order accurate preconditioned AUSM scheme is used for Cartesian grids whereas Roe's scheme is applied on nearbody grids. An implicit time integration using self-adaptive Newton method [START_REF] Daude | Self-Adaptive Newton-based iteration strategy for the LES of turbulent multi-scale flows[END_REF] is applied, with a maximum of 20 subiterations per time step. Currently, rotor blades are resolved by nearbody structured grids in relative motion with the Cartesian background grids. Unsteady Chimera transfers are applied at each subiteration of the time integration.

The fuselage can be either represented by a structured body-fitted mesh overlapping the off-body Cartesian mesh or by an immersed boundary within the Cartesian mesh. As the objective of this paper is to evaluate and compare Immersed Boundary Methods developed in Helios/ROAM and FAST/IBM for the capture of interactional aerodynamics of a rotorcraft configuration, the fuselage is represented by an Immersed Boundary.

FAST/IBM approach is a Ghost Cell Direct Forcing formulation, designed for attached turbulent boundary layers, as de-scribed in details in P éron et al. [START_REF] P Éron | An immersed boundary method on Cartesian adaptive grids for the simulation of compressible flows around arbitrary geometries[END_REF] . This approach consists in marking cells as fluid or solid and forcing the solution at the fringe of solid points, in the close vicinity of the obstacles, such that the velocity is implicitly equal to zero at an immersed boundary describing a no slip boundary condition for instance.

In figure 2, the immersed boundary (black curve) tiles the domain into a solid region (left-hand side of the black curve) and fluid region (right-hand side). IBM forced or target points are the two layers of cells at the fringe of solid points (red dots). The IBM reconstruction needs an IB wall point (that is the projection of an IB target point onto the surface describing the immersed obstacle) and an IB image point (blue dot), the three of them being aligned along the wall normals. To update the solution at IBM target point A in red at a given iteration n, the solution at point B is used. As this point is not a discretization point, its solution is interpolated from its neighboring computed cells (donor points as green dots). Then, pressure and density defined at IB im-age point B are set at IB target point. In the case of a turbulent boundary layer, the velocity at point A is obtained as follows: the friction velocity is computed at IB image point B and set at point A consequently. Then Musker's wall function is applied to reconstruct the tangential velocity at IB target point A. This wall function is algebraic and fits both the log layer and the viscous sublayer. The normal velocity component is obtained by a linear reconstruction directly from the normal velocity at point B. The Spalart-Allmaras pseudo-viscosity at IB target points is also reconstructed (details can be found in the appendix of the paper 6 ) .

4.SIMULATION OF THE ROTOR-FUSELAGE INTERACTION OF THE D365N CONFIGURATION

This section compares mid-fidelity simulations using FAST/IBM and Helios/ROAM with experimental results and three high-fidelity simulations. The three high-fidelity results are from Potsdam using OVERFLOW and Renaud using elsA in 2009 [START_REF] Potsdam | Unsteady Computations of Rotor-Fuselage Interaction[END_REF] , using overset, structured curvilinear grids with low-Mach preconditioning for the fuselage. The last high-fidelity result is a recent result obtained using Helios mStrand without low-Mach preconditioning. The FAST/IBM simulation performed by ONERA is a URANS simulation using a wall-modeled IBM onto the fuselage and overset body-fitted grids for the blades. U.S. Army approach relies on Helios/ROAM, where the fuselage is also modeled by an immersed boundary, but without any wall function and an actuator line model for the rotor. A same spacing has been chosen for both simulations near the immersed boundary within the rotor wake in order to be comparable.

4.1.Grid description

ROAM and FAST/IBM simulations rely on the same spacings in the vicinity of the fuselage and around the rotor disk in order to compare the IBM approaches with the same near-wall resolution. A close-up view of the skeleton of the Cartesian mesh is represented in Figs. 3(a Helios/ROAM skeleton tree is made of 33,386 octants, whereas FAST/IBM skeleton is made if 45,531 octants; indeed, the refinement levels close to the fuselage must be wide enough in the normal direction to the obstacle to maintain forced and image points in the same refinement level. Each octant is then filled with 16 points per direction, resulting in final Cartesian meshes made of 136.7 and 186.5 million points respectively. The finest cell is 0.977 m in the vicinity of the hub and the tail to represent sufficiently the sharp features of the geometry by the IBM. In other regions close to the fuselage, the cell spacing is twice as big. For FAST/IBM approach, the actual distance to wall in wall unit y `is around 100. For FAST/IBM simulation, the four main rotor blades are meshed by a set of structured grids of 2 million points, overlapping the Cartesian mesh (Fig. 4 ).

Figure 5 shows the actuator line representation of the rotor blades results in clear high-and low-pressure regions in the background mesh, as shown by the red and blue areas.

4.2.Rotor loads and wake visualization

The rotor thrust for Helios/ROAM and FAST/IBM simulations is compared to highfidelity results using Helios with mStrand (body-fitted solver) and reference overset results from Potsdam et al in 2009 [START_REF] Potsdam | Unsteady Computations of Rotor-Fuselage Interaction[END_REF] . Experimental results from the F1 tunnel are also included for comparison [START_REF] Potsdam | Unsteady Computations of Rotor-Fuselage Interaction[END_REF] .

Time-averaged rotor thrust values are summarized in Table 1. FAST/IBM was run using two different time steps, the first one corresponding to a variation of azimuth per time step ∆Ψ " 0.5 ˝and the second one to ∆Ψ " 0.25 ˝. It has been observed that the variation of azimuth ∆Ψ " 0.5 ˝leads to an underprediction of the thrust. Instead using ∆Ψ " 0.25 ˝per time step seems to be a good compromise regarding the accuracy, which was also observed by Potsdam et Figure 6 shows the rotor thrust over a full revolution of the rotor with averages shown as dotted lines. Despite a good agreement between all the simulations, they underpredict the experimental thrust of 10% roughly. The discrepancy may be due to assumptions made in HOST when trimming the rotor to obtain deflections for CFD. For consistency with previous high-fidelity results from 2009, the same rigid blade deflections were used and the rotor was not re-trimmed using CFD-CSD coupling. For the purposes of this work comparing mid-and high-fidelity results is more relevant than matching experimental results.

Figs. 7(a) and (b) display an isosurface of the Q criterion at Q " 11, 548 s ´2, highlighting the impingement of the rotor wake on the fuselage. Compared to FAST/IBM, which has boundary-layer-resolved meshes, the actuator line method in Helios ROAM results in larger, diffused vortices trailing from the blade tips. The IBM methodology in Helios ROAM without a wall model results in flow separation on the fuselage and many small-scale flow features hitting the tail region. These small flow features are likely not physical, since they do not appear in the results from FAST/IBM or the high-fidelity results from Helios/mStrand or Potsdam et al..

Another deficiency of the ROAM actuator line method is poor resolution of the nearwake region behind the blade. At the same Q-criterion level the body-fitted blades solution from FAST/IBM shows a vortex sheet behind the blade. Since ROAM uses lifting line theory at discrete control points along the blade, the injection of source terms in the background mesh is not adequate to resolve span-wise gradients that result in the nearwake trailed vorticity.

The rotor and fuselage wake can be studied in further detail by analyzing vorticity at different planes of the background domain. Two planes were previously used for comparison against experimental PIV data in the work of Potsdam et al.. Vorticity in the first plane, located at X " 1.52m behind the aircraft nose or (1.0125m behind the hub center), is shown in Fig. 8 for FAST/IBM and ROAM. The vortex structure is clearly asymmetric, as the blade root vortices visible in Fig. 7(a) coalesce on the left side of the fuselage and impact the horizontal tail. The ROAM solution from Fig. 8(b) shows many small-scale vortex structures behind the tail, as previously noted. The large rotor wake vortices on both the advancing and retreating sides of the aircraft are well captured by both solvers. The advancing (left) side has a stronger nearwake vorticity since the loads are distributed across the entire blade span. On the retreating (right) side the loads are more concentrated near the tip as a result of low relative velocity on the inboard sections of the blade. Since FAST/IBM has fully resolved blade meshes, the wake solution contains more clearly defined vorticity.

Figure 9 shows a Y " 0.475 slice near the blade mid-span with contours of Y ´vorticity on the retreating side of the aircraft. Again vorticity from Helios ROAM appears larger and more diffused compared to the tighter vortex structures from FAST/IBM. As the vortices convect behind the aircraft, the structures in FAST/IBM dissipate whereas with Helios/ROAM they maintain their strength longer. This is likely a result of the 2 nd order spatial scheme used in FAST/IBM compared to the 5 th order convective scheme in Helios.

5.CONCLUSIONS

In this paper, an immersed boundary method (IBM) is assessed for the interactional aerodynamics between the rotor and the fuselage of a Dauphin helicopter in the frame of the US/France agreement on helicopters AHFIR.

This method has been developed recently both by ONERA and US Army on adaptive Cartesian grids, within FAST/IBM and Helios/ROAM frameworks respectively. This method enables to overcome the mesh generation effort as the mesh does not need to conform to obstacles; an immersed boundary condition (IBC) is applied to reconstruct implicitly the wall boundary condition. CFD simulations are performed using a Cartesian solver that has proven efficient in terms of memory requirements and CPU time in comparison with curvilinear solvers.

Immersed boundary approaches are a good compromise between standalone lowfidelity tools and high-fidelity methods, the first ones usually lack the resolution necessary to capture interactions with the fuselage. Coupling low-fidelity tools with CFD for the fuselage to improve the accuracy of the prediction of the download can be achieved at the expense of generating a body-fitted mesh, which could be a tedious task for new configurations.

For that purpose, US Army approach consists in combining an actuator line model for the rotor with an immersed boundary for the fuselage within Helios/ROAM, whereas ONERA approach relies on a wall-modeled immersed boundary method coupled with overset grids for the rotor blades. It has been demonstrated that the interactional effects can be captured using both approaches: both approaches demonstrate the ability to capture the skin pressure coefficient over different parts of the body, showing reasonably good agreement with both the experiments and more computationally expensive wall-resolved body-fitted CFD results. Analysis of the wake, especially in the vicinity of the empennage, demonstrates the capture of the interactional aerodynamics of this configuration. Helios/ROAM is not predicting the 4-per-rev signature of the fuselage drag as well as FAST/IBM. This is probably due to the immersed boundary method currently used ROAM which does not model the turbulent boundary layer, as opposed to the FAST/IBM approach which does through the application of a wall function at the IBM points on the fuselage.

In order to have at disposal tools of different fidelity, ONERA aims also at introducing a low-fidelity tool such as the actuator line model in ROAM to be combined with IBM. Future work will consist in applying the immersed boundary method onto a main rotor hub, for which the method should be wellsuited as this kind of geometry can be complex and thus tedious to mesh using bodyfitted approaches. The objective is to assess the capability of the IBM to predict the turbulent separated flows around geometrically complex parts of the rotorcraft. These validation studies will open the path to an efficient tool for helicopter download prediction of new configurations, with a good compromise between accuracy and process time.

  Gaussian injection of source termsin background mesh equations.

Figure 1 :

 1 Figure 1: Illustration of actuator line implementation to simulate a rotor blade using the injections of source terms from sectional loads.

Figure 2 :

 2 Figure 2: FAST/IBM approach: type of points used in the IBM reconstruction.

Figure 3 :

 3 Figure 3: Octree skeleton of the Cartesian mesh for (a) FAST/IBM and (b) Helios/ROAM

Figure 4 :

 4 Figure 4: View of FAST/IBM mesh for the D365N configuration: four structured grids for each blade with the Cartesian mesh.

Figure 5 :

 5 Figure 5: Helios/ROAM pressure contours (skeleton of the Cartesian mesh is represented).

Figure 6 :

 6 Figure 6: Comparison of rotor thrust between Helios/ROAM and FAST/IBM calculations against high-fidelity simulations and experiments.

Figure 7 :

 7 Figure 7: Isosurface of Q-criterion Q " 11548 s ´2 and pressure contours for (a) FAST/IBM and (b) Helios/ROAM.

Figure 8 :

 8 Figure 8: CFD wake X´vorticity (s ´1) on plane X " 1.52m for (a) FAST/IBM and (b) Helios/ROAM.

Figures 10 (

 10 Figures 10(a) and (b) show the pressure reconstructed on the immersed surface and averaged over a quarter-revolution. The pressure is interpolated from the background mesh for ROAM, and using a Moving Least Squares reconstruction of pressure values at IBM target points for FAST/IBM. The centerline pressure coefficient is extracted and plotted along the upper-surface of the aircraft in Fig. 11. Stagnation pressures appear as red contours in the contour plot or spikes to C P « 1 in the centerline plot. Stagnation points appear in expected locations on the leading-edge surfaces of the fuselage. At the vertical tail, however, the pressure coefficient increases to well above 1, indicating an interaction from the main rotor.

  Figure 10: Time-averaged pressure coefficient on the fuselage.

Figure 11 :

 11 Figure 11: Time-averaged centerline pressure coefficient.

Figure 12 :

 12 Figure 12: Comparison of fuselage lift between different approaches.

Figure 13 :

 13 Figure 13: Comparison of fuselage drag between different approaches.

  Helios are implemented in a module called the Reduced Order Aerodynamic Model (ROAM). ROAM is used in Helios alongside a Cartesian solver, either SAMCart or Orchard. For this work, only the octree-based AMR solver Orchard 5 is used as the background Cartesian solver.

Table 1 :

 1 al. for the high-fidelity simulations. All Helios simulations used a timestep of ∆Ψ " 0.25 Comparison of average rotor thrusts. The first four approaches are from the current work; the last three are from referenced work.

			˝.
		∆Ψ	Z b
	Helios ROAM	0.25 ˝13.53
	Helios mStrand 0.25 ˝13.48
	FAST/IBM	0.5 ˝13.2
	FAST/IBM	0.25 ˝13.44
	F1 exp	-	14.5
	AFDD 2009	0.25 ˝13.5
	ONERA 2009	0.25 ˝13.36

6.ACKNOWLEDGMENT

The authors would like to acknowledge the support of the France-USA AHFIR cooperation represented by Arnaud Le Pape from ONERA side and Mark Potsdam and Rohit Jain from US Army side.


D Zint

N Maruani

P Alliez

M Rouxel-Labbé

M Rouxel

Feature-Preserving Offset Mesh Generation from Topology-Adapted Octrees

Keywords:

Input shape (b) Dual Contouring (2%) (c) Offset (2%) (d) Dual Contouring (5%) (e) Offset (5%) Figure 1: Method overview: (a) The input consists of a mesh and two main user-defined parameters: offset radius and maximum normal deviation; (b,d) The topology of the initial offset mesh is computed with Dual Contouring on a topology-adapted octree; (c,e) A high-quality offset mesh is obtained through remeshing.

Introduction

Rolling a ball of fixed radius over a surface describes a volume whose boundary is referred to as the offset surface. More formally, this operation is known as the Minkowski sum of the surface and the ball. The offset surface can also be defined as the level-set of the unsigned distance field to the input. The construction of offset surfaces is a fundamental tool in many applications such as computeraided design (CAD), collision detection, path planning, boundary layer mesh generation, architectural design, design exploration, etc. Despite being a simple, well-defined mathematical operation, the reliable generation of discrete offset surfaces is a notoriously difficult scientific challenge, as even simple inputs can lead to complex offset surfaces (see Figure 2). This complexity has caused offset surface computation to be studied mostly in the context of discrete inputs and outputs. Using discrete inputs and outputs is nevertheless reasonable, given the discrete nature of most real-world data. In addition, a discrete approximation of the offset that lies within a user-defined error bound is often sufficient, or even required.

Focus and Problem Statement

This paper addresses the problem of generating a discrete offset in the form of a surface triangle mesh, given a surface mesh as input, with two main user-defined parameters: an offset radius δ and a maximum normal deviation in degrees σmax. Minimum edge length and maximum octree depth are optional parameters that provide additional control to the user. The maximum normal deviation σmax provides a means to control the approximation error. The minimum edge length and the maximum octree depth limit the mesh complexity. Our main focus is to guarantee a closed, combinatorially 2-manifold output surface triangle mesh through a process that is robust to input defects such as self-intersections or holes. We also aim to generate offset meshes with well-shaped (isotropic) triangles everywhere, except near sharp features subtending small angles.

Related Work

The literature abounds with the generation of offset meshes as it is a complex, multifaceted problem. We observe no consensus on the best methodological approach due to multiple dilemmas: exact vs. inexact, over-refinement followed by filtering vs. coarse-to-fine, etc. As such, some approaches aim at constructing data structures in relation to the structure of the distance function in space (generalized 3D Voronoi diagrams), while others aim at understanding the structure of the exact local offset. Some approaches reformulate the problem as the approximation of Minkowski sums, and many approaches utilize volumetric discretizations.

Generalized 3D Voronoi Diagram. The distance field d(x) of a triangle mesh is C∞ continuous almost everywhere but is only C 0 in some areas. More specifically, the distance field is only C 0 whenever the nearest primitive (vertex, edge, face) changes. For example, when we move in the proximity of a concave edge, the gradient of the distance field jumps at some point. The nearest primitive to each point in space is encoded by the generalized 3D Voronoi diagram. The distance field is C∞ inside Voronoi cells, while it is only C 0 on its bisectors. Offset generation can thus be seen as a sub-problem of generating the generalized 3D Voronoi diagram. If we had the exact 3D diagram, we could simply compute the offset for each cell, a simple task due to the C∞ continuity, and then cut the offsets at the bisectors. However, contrary to 2D where generalized Voronoi diagrams are well studied, with available implementations [START_REF] Karavelas | 2D Segment Delaunay Graphs[END_REF] . This elegant approach makes it possible to compute the exact topology of the offset around a vertex (with the crucial requirement of computing the exact generalized Voronoi diagram on the sphere) and therefore to create an offset mesh with arbitrary precision. However, this method does not yet extend to global intersections, i.e. when two non-adjacent elements are closer than twice the offset distance.

Minkowski sum. Varadhan and Manocha approximate the Minkowski sum of polyhedral models by computing distance fields and reconstructing the surface with a variant of Marching Cubes [START_REF] Varadhan | Accurate Minkowski sum approximation of polyhedral models[END_REF]. Campen and Kobbelt proposed a powerful method to compute intersections of polygonal meshes and leverage it to generate approximate Minkowski sums between the input mesh and a sphere [START_REF] Campen | Exact and robust (self-) intersections for polygonal meshes[END_REF][START_REF] Campen | Polygonal boundary evaluation of Minkowski sums and swept volumes[END_REF]. A polygonal surface mesh approximates the sphere swept over the input mesh. For the context of manufacturing, approximating this sphere by a zonotope is often sufficient. Martinez et al. use this property to provide an efficient method for approximating offset surfaces [START_REF] Martinez | Chained segment offsetting for ray-based solid representations[END_REF]. However, when aiming for high accuracy, especially in concave areas, the sphere must be discretized with a dense mesh.

Volumetric discretizations. To avoid the difficulty of constructing exact offset surfaces, most methods shift their objective towards submitted to Eurographics Symposium on Geometry Processing (2023) approximating the offset. Marching Cubes [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF] and Dual Contouring [START_REF] Ju | Dual contouring of hermite data[END_REF] are both widely applicable isosurfacing methods that have been constantly improved over the past decades. Qu et al. use an approach similar to Schaefer and Warren's Dual Marching Cubes that is based on a structured irregular grid [QZS *04]. While this approach can represent sharp features, it is based on a dense voxel grid, which limits its practicality. Liu and Wang successfully use a modified version of Dual Contouring to generate intersectionfree offset surfaces [START_REF] Liu | Fast intersectionfree offset surface generation from freeform models with triangular meshes[END_REF]. This method is also based on a dense voxel grid. Pavić and Kobbelt use an octree to avoid the overhead caused by voxel grids [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF]. However, the octree is subdivided to the lowest level whenever cells are non-empty, causing milder but similar types of over-refinement. Features are reconstructed in a post-processing step by setting a threshold on the normals and adding features by subdividing faces and edges. Outliers and lowquality elements are removed by further post-processing. Although the aforementioned methods are reliable and robust, the uniform grid resolution (or maximal octree depth) aspect hampers scalability and generates transient over-refined meshes.

GPU-based and sampling methods. Whenever a method relies on structured data, GPUs can reduce runtime significantly by performing several thousands of operations in parallel. Wang and Manocha use Layered Depth Images on GPUs to sample the surface of the input mesh [START_REF] Wang | GPU-based offset surface computation using point samples[END_REF]. This method can deliver more than 100 times speedup in comparison to related CPU methods. Chen et al. compute offset surfaces with an efficient dexel data structure that stores an array of cells containing a list of intersections with rays [START_REF] Chen | Half-space power diagrams and discrete surface offsets[END_REF]. While these methods solve the runtime issues that many volumetric approaches have, they still generate meshes with large complexity. Meng et al. distribute sample points uniformly across the offset [MCS*18]. This method requires dense sampling to ensure the correct representation of concave creases.

Meshing and Remeshing. In many cases, the raw output of the methods described above requires post-processing to treat issues like self-intersections or low-quality elements. Additionally, remeshing can be used to recover features and reduce the discretization error. One such technique is the isotropic remeshing introduced by Botsch and Kobbelt, which focuses on element quality [START_REF] Botsch | A remeshing approach to multiresolution modeling[END_REF]. While this method succeeds in its objective, it is not sufficient in our context as the method assumes that the underlying surface is continuously differentiable, i.e. it does not contain sharp features, which is not true for offsets. Shen et al. approximate implicit surfaces in [START_REF] Shen | Interpolating and approximating implicit surfaces from polygon soup[END_REF] but their method also smoothes away small features. The remeshing approach contributed by Pavić and Kobbelt in [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF] first detects feature lines in the previously generated offset mesh by setting a threshold on the normal deviation within a triangle and reconstructs them by adding vertices along the feature line and flipping edges to reconstruct the features. Smoothing is applied to improve element quality and reduce the discretization error.

Method Overview

Given an input triangle surface mesh T and an offset radius δ, the offset surface S(T , δ) is defined as the isosurface d(x, T ) = δ where d(x, T ) is the distance field of T . If the input is a closed mesh, the distance field of T can be signed, and the offset is composed of multiple closed surfaces. We denote by V(T , δ) the volume bounded by S(T , δ).

Our approach consists of two distinct yet complementary steps, depicted in Figure 1:

Octree construction and Dual Contouring: This step is designed to generate an initial offset mesh with the desired topology: closed, combinatorially 2-manifold, and approximately (up to a user-defined tolerance) equivalent to the offset surface. The octree is constructed using several topological criteria to drive refinement. The geometry of the generated mesh is still crude with, e.g., no particular care given to sharp features. See Section 4.

Remeshing: This step applies a novel remeshing algorithm tailored to offset surfaces to improve the geometric precision while keeping the topology unchanged. See Section 5.

Positioning and Contributions

Our approach builds upon several existing methods related to offset mesh generation and sharp feature recovery.

Adapted Octrees. At first glance, the method of Pavić and Kobbelt [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF] is closest to ours: it is based on a similar twophase approach. Its first phase is also based on an octree, but the octree is refined to its maximum depth everywhere around the offset surface. This has a number of disadvantages. First, it causes a significant computational overhead because one subdivides cells even in topologically and geometrically simple regions, see Section 6.2.1. It also means that the maximum depth is limited as the octree becomes exceedingly large. Furthermore, the user must select a maximum octree depth. While the relationship to the approximation error is clear, the one to the topological errors is not intuitive. Finally, such a dense refinement translates into a complex intermediate mesh, which slows down further processing steps such as remeshing.

Varadhan et al. make use of a topology-adapted octree for computing Minkowski sums [START_REF] Varadhan | Accurate Minkowski sum approximation of polyhedral models[END_REF][START_REF] Varadhan | Topology preserving surface extraction using adaptive subdivision[END_REF]. The offset needs to be approximated by sweeping a discretized sphere over the input mesh. The offset quality depends on the sphere discretization and the computation of pairwise-convex Minkowski sums is computationally expensive.

Our first main contribution is an octree generation that estimates the offset conservatively at all times and is therefore guaranteed to never miss the exact offset, up to floating point precision. Additionally, our conservative cell-to-triangle distance approximation is computationally lighter than the methods presented in [START_REF] Varadhan | Accurate Minkowski sum approximation of polyhedral models[END_REF] and [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF], and does not require the input to be closed or manifold, e.g. see Figure 3. Our topology-adapted octree reduces computational overhead and output complexity by performing fewer subdivisions in regions with disk topology. In our framework, the maximum depth parameter is not as much a limiting factor as in other methods because the maximum depth can be set to a large value without hampering runtime when such a maximum depth is not utilized. The intermediate offset mesh is reconstructed by Dual Contouring, which generates closed meshes by design. By post-processing this output, we also guarantee that the final offset mesh is not just closed but also topologically 2-manifold everywhere.

Remeshing. Meshes generated by Dual Contouring may contain low-quality or tangled elements. Our second main contribution reduces this limitation and improves on existing remeshing algorithms. More specifically, we contribute a novel remeshing algorithm tailored to offset surfaces, which generates elements with overall good quality while ensuring that the offset mesh still approximates the offset well. The remeshing algorithm departs from earlier approaches in two ways: (1) It is driven by normal deviation instead of by the common edge length or approximation error, and (2) Its vertex relocation operator utilizes the quadric error metric (QEM) initially devised for mesh decimation [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]. The resulting remeshing algorithm is feature-aware and recovers a faithful discretization of concave creases and corners. It also yields highquality reconstructions of challenging geometric features such as thin tunnels. Finally, it is independent of the octree and can also be applied as a post-processing step to other offset mesh generation methods such as [PRH*22] or [START_REF] Campen | Polygonal boundary evaluation of Minkowski sums and swept volumes[END_REF].

QEM-based remeshing has been explored before, for example by Valette et al. [START_REF] Valette | Generic remeshing of 3D triangular meshes with metricdependent discrete Voronoi diagrams[END_REF]. This is however a fundamentally different remeshing approach: surfaces are sampled and sampling points are regrouped into clusters. Starting from those clusters, a Centroidal Voronoi Diagram is generated by minimizing an energy term. Error quadrics are used to relocate centers of Voronoi cells, which has the effect of snapping vertices to feature lines or corners. Although effective, this method cannot be applied to offset surfaces due to its implicit nature. If we were to approximate the offset with a mesh, we could not guarantee that this mesh does not contain any tangling that cannot be resolved by remeshing. In addition, the resulting mesh might contain non-manifold edges or vertices. Finally, in the case of Valette et al., the method generates a user-defined number of vertices. While this is an important feature for surface decimation, when generating offset meshes the complexity of the output is often unknown a priori.

Octree Construction and Dual Contouring

The octree construction begins with its initialization by a single cell made large enough to contain the input mesh and its offset. We utilize the offset surface S(t i , δ) and the enclosed volume V(t i , δ) of individual triangles of T = {t i }. Each octree cell is subdivided recursively until one of the following criteria is met:

1. No triangle offset V(t i , δ) intersects the cell; 2. The cell is fully enclosed inside the offset surface of a triangle V(t i , δ); 3. The offset surface S(T , δ) is topologically equivalent to a disk within the cell; 4. A user-defined maximum octree depth is reached.

The first two criteria are devised to stop the refinement of cells that do not intersect the offset. They are regrouped in the so-called Intersection criterion, described in Section 4.1. The third criterion aims to stop the refinement of cells that are intersected by the offset, but where subdivision does not add any relevant topological information. We call it the Disk criterion, Section 4.2. Finally, we might have to subdivide cells to circumvent the known issue that Dual Contouring may generate non-manifold vertices and edges. This step is detailed in Section 4.3. All criteria are evaluated conservatively, meaning that we might over-refine but never terminate too early.

Checking these criteria with all triangles at all cells would be prohibitive. Instead, we store for each octree cell all triangles such that V(t i , δ) intersects the circumscribing sphere of the cell, similarly to Pavić and Kobbelt [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF]. Our initial cell thus contains all input triangles and each new subdividing cell needs only find those that are relevant within the triangles of its parent cell.

Intersection Criterion

Octree cells are only of interest if the offset surface intersects them. However, checking for the exact intersection of a cell with all triangle offsets is costly. We accelerate this test by approximating the cell with its enclosing sphere. To check whether the offset of a triangle t i intersects this sphere, we find the nearest point pi submitted to Eurographics Symposium on Geometry Processing (2023) of the cell center pc on the triangle. Next, we compute the distance d i = pipc to that point. The sphere intersects the triangle offset if the difference between the absolute offset radius |δ| and the distance is at most the radius r of the sphere,

d i -r < |δ| < d i + r. (1)
Triangles whose offset do not intersect the sphere are culled from the cell's triangle list.

As cells are divided, some cells become fully enclosed within the offset of a triangle,

d i + r < |δ|. (2)
In this case, the whole triangle list is discarded and the cell is removed from the subdivision queue.

So far, the intersection criterion only considers the unsigned offset radius. When the input mesh is closed, we use the signed offset by checking whether the cell contains the correct offset sign for any of its primitives,

d s,i -r < δ < d s,i + r. (3)
We utilize the signed distance d s,i with a negative sign when the cell center pc is on the bounded side of a closed mesh. The sign is determined with Side_of_triangle_mesh from the CGAL component Polygon Mesh Processing [START_REF] Loriot | Polygon Mesh Processing[END_REF]. A negative offset radius corresponds to offsetting towards the inside, see Figure 4. A cell in which no triangle t i satisfies Equation (3) is not subdivided any further, as it does not contain the offset of interest.

Disk Criterion

An offset tends to have simple topology in most areas but it can be very complex in areas where multiple offset sheets touch. Capturing the topology then requires a high level of octree refinement.

As we wish to preserve the topology of the exact offset surface, we apply a criterion to capture topological features: A cell in which the offset surface is topologically equivalent to a disk is no longer subdivided. To evaluate this predicate, we construct a sphere with radius δ for every projection point pi from Section 4.1. If the set of spheres has a non-empty intersection, the union of triangle offset surfaces ti∈T V(t i , δ) form a star domain, as the offset of a triangle is always convex, see Figure 5. Any star domain is simply connected and therefore of genus 0. It follows that any closed and manifold subset has disk topology.

A non-empty intersection exists if there is a point with a distance to all primitives smaller than the offset. Such a point is found by constructing the minimal sphere that contains all projection points. If this sphere has a radius smaller than the offset, then an intersection exists. We also perform a cheap test by checking if all triangles share a vertex, which trivially guarantees the existence of a non-empty intersection.

We need to ensure that the above genus 0 object can also be reconstructed with Dual Contouring. Therefore, if all cell corners are located on the same side of the offset, the cell must be further subdivided.

While we can guarantee that the union of triangle offset surfaces forms a star domain, the second part, showing that the intersection of the cell with the star domain is a closed manifold, is more complex. A proof that the intersection of all triangle offset surfaces and also the cell is non-empty would be required. In the case of creases subtending very small angles, Dual Contouring tends to yield the wrong topology, as the second part of the disk criterion is no longer met. We alleviate this issue by subdividing a cell if the offset normals within the cell vary by more than 120 degrees. Nevertheless, when the offset contains very sharp creases, we cannot guarantee the topological correctness of our offset mesh as it may contain small holes or disconnected islands.

Our disk criterion relates to the subdivision criteria introduced by Varadhan and Manocha [START_REF] Varadhan | Topology preserving surface extraction using adaptive subdivision[END_REF], but its use and evaluation differ significantly.

Subdivision for Manifoldness

A known pitfall of the Dual Contouring method is the possible creation of non-manifold edges and vertices. If an occurrence of such a configuration is detected, we apply a subdivision step to the cells involved. If such a subdivision is forbidden because we reach the user-defined maximum octree depth, non-manifoldness is resolved by duplicating vertices that are either non-manifold or incident to an edge that is non-manifold. Alternatively, replacing Dual Contouring by the Manifold Dual Contouring presented in [START_REF] Schaefer | Manifold dual contouring[END_REF] would also be feasible.

Offset-Aware Remeshing

We use the common mesh operators used for surface remeshing: edge splits, edge collapses, edge flips, and vertex relocation. Instead of driving the remeshing process by the edge length as common in previous work, we utilize the normal deviation between the offset mesh triangles and the exact offset.

Definition 5.1 The normal deviation σ(t) of a triangle t is the maximal angle between the triangle normal nt and the offset normal no in the offset area this triangle is mapped to:

σ(t) = max ∠(nt , no(x)),
where x denotes any location in the triangle offset area.

Normal deviation can be understood as a measure of offset curvature, but with the advantage that normal deviation is still well defined at concave creases, where the offset curvature is discontinuous. Therefore, it can be used to isotropically refine smooth regions, implicitly adapting elements to the curvature, and detect concave creases.

Each operation fulfills a different task in the remeshing. Edge splits increase the number of vertices in regions where the normal deviation is too large. More vertices lead to more flexibility in adapting the mesh to the desired geometry. Edge collapses do the exact opposite, they reduce the mesh complexity in regions with low normal deviation, e.g. in planar regions. Edge flips improve element quality and snap edge to concave creases. Vertex relocation minimizes the normal deviation using QEM. Operations are performed subsequently, in the order they are mentioned above.

Computing Normal Deviation

Computing the exact normal deviation is time-consuming as it requires the use of optimization for searching the local mapping. As the distance field is only C 0 -continuous, the optimization problem is challenging. Instead, we sample the triangles of the offset mesh uniformly and project the sample points onto the input mesh. For each sample point, the offset normal nt is then constructed by the normalized vector pointing from the projection point ps to the sample point ps. The maximum normal deviation of a triangle is computed by taking the maximum angle between the triangle normal and all offset normals at the sample points within the triangle. We use CGAL's AABB tree data structure [START_REF] Alliez | 3D Fast Intersection and Distance Computation[END_REF] to project sample points onto the input mesh.

Edge Split

We perform a batch of edge split operations. The list of candidate edges for splitting is found by searching for triangles with σ(t) > σmax, where σmax is a user defined maximum normal deviation. The longest edge of such a triangle is considered as a candidate. The new vertex inserted by the edge split operator is located as described in the vertex relocation step, see Section 5.5.

Note that we perform only a single batch of edge split operators, without considering the newly created edges as candidates for recursive splitting. The edge split operator does not guarantee a reduction of normal deviation, but adding more vertices to the mesh adds new degrees of freedom that are later used in the vertex relocation.

Our method contains an optional parameter so that the user can set a minimum length. If an edge is already smaller than the minimal length, the edge will not be refined. This option avoids resolving smaller features that are not of interest to the user.

Edge Collapse

In areas that contain more vertices than necessary, we perform a series of halfedge collapse operations to reduce the mesh complexity. A vertex is collapsed into another one when the normal deviation in all its incident faces is smaller than the user-defined maximum deviation.

Our remeshing operations were designed under the assumption that the mesh is isotropic. Therefore, we must avoid rapid changes in element size. If we collapse all edges only according to the normal deviation, flat areas would be fully collapsed, and the assumption of isotropy would not hold any longer. We avoid this by only collapsing edges that are smaller than twice the maximum edge length in convex offset regions, lmax = 2δ sin(σmax).

Edge Flip

Edge flip operators are used to favor well-shaped (isotropic) triangles. As they should not interfere with the optimization for normal deviation, we introduce the following metric for a triangle:

mq = 2 √ 3A/(l 2 1 + l 2 2 + l 2 3 ) (4) m ∠ = 1 -∠(nt , no)/90 • (5) m flip = mqm 3 ∠ , (6)
where A denotes the triangle area, l i , i ∈ [1, 2, 3] is the length of each edge incident to the triangle, nt is the triangle normal, and no is the offset normal at the triangle center. Equation (4), also known as the mean-ratio metric, was already used for mesh optimization [START_REF] Bank | Pltmg: A software package for solving elliptic partial differential Equations: Users[END_REF][START_REF] Rangarajan | Provably robust directional vertex relaxation for geometric mesh optimization[END_REF]. Equation (5) denotes the normalized angle between the triangle and the offset surface normal at the triangle center. Equation (6) combines these two metrics, with a greater preference given to the angle metric. Finally, an edge is flipped if the minimum m flip quality of its two incident triangles is improved. Flips on triangles with σ(t) < σmax that causes the normal deviation to exceed the user-defined maximum are prohibited.

Vertex Relocation

Vertices are relocated using a combination of Laplace smoothing and minimization of error quadrics. This approach is similar to the one presented in [START_REF] Valette | Generic remeshing of 3D triangular meshes with metricdependent discrete Voronoi diagrams[END_REF]. Intuitively, our vertex relocation operator improves the shape of triangles in smooth areas and snaps vertices onto creases or corners if there are any within the umbrella of a vertex.

Given a vertex v at position x, we first compute the center of mass of its neighboring vertices,

x = 1 |N(v)| ∑ xn∈N(v) xn, (7)
submitted to Eurographics Symposium on Geometry Processing (2023)

where N(v) denotes the set of neighboring vertex locations and |N(v)| denotes the number of neighbors. This step is referred to as Laplace smoothing [START_REF] Field | Laplacian smoothing and Delaunay triangulations[END_REF].

We then compute the weighted error quadrics of tangential planes at all sampling point projections onto the offset ŝ. We use the same sampling as for computing the normal deviation, see Section 5.1, on each triangle incident to v. The tangential planes are weighted by the area of the triangle they belong to. For each point ŝi , the offset normal ni and triangle area a i are required. The sum of weighted error quadrics is minimized by solving the following linear system of equations:

Ax = b (8) A = n samples ∑ i=1 a i • ni n i (9) b = n samples ∑ i=1 a i • n i ŝi . ( 10
)
Solving Equation ( 8) is often not possible because A does not have full rank, e.g. in planar regions all normals ni are similar. We compute the pseudo-inverse with singular value decomposition, as proposed by Lindstrom [Lin00]:

A = UΣV (11)
A + = VΣ + U (12) x = x + A + (b -Ax ). (13)
Using the pseudo-inverse requires a starting position x , chosen as the previously computed center of mass, Equation (7). This way, the vertices are relocated by the error quadrics minimization only when there is a distinct minimum. If quadrics minimization yields a line or plane of optimal points, then we choose the one that is closest to x .

Results

We evaluate our method on a wide range of input meshes with varying offset radii. Section 6.1 demonstrates robustness by running our method on the Thingi10k dataset. Section 6.2 offers a more in-depth analysis, discussing the impact of parameters on quality and runtime.

Our method was implemented in C++ and depends on the libraries Eigen [GJ*10] and CGAL [START_REF][END_REF]. For measuring runtime, the code was serially executed on a cluster with Intel Xeon Platinum 8268 processors that have a base frequency of 2933 MHz.

Validation

We generated offsets for all meshes from the Thingi10k dataset [START_REF] Zhou | Thingi10K: A Dataset of 10,000 3D-Printing Models[END_REF]. Unless stated differently, all results presented in this paper use the following default values. We use a maximum octree depth of 10 plus two additional levels for resolving non-manifoldness. The remeshing step performs 10 iterations and uses a maximum normal deviation of σmax = 7 • . Figure 6 showcases results generated from the Thingi10k dataset.

Our implementation expects non-degenerate triangles as input, but the underlying principles make it feasible to operate on point clouds or meshes with dangling edges. Eventually, we ran our method on 9,952 meshes. All generated offset meshes are closed manifolds. Besides the aforementioned default values, we set a relative offset radius of 10% and 5% with respect to the largest edge of the bounding box, see Figure 7. For evaluation, we densely sample the offset mesh and compute the distance to the offset and the normal deviation using the Libigl geometry processing library [JP*18].

We denote by εm the discretization error, measured as the mean distance of a mesh to the offset surface relative to the given offset radius, and denote by σm the mean normal deviation. For almost all offset meshes, εm is below 0.5%. This means that for a mesh with an offset radius of 0.05, the exact offset is missed by an average of 0.00025. The average maximum distance of all offset meshes relative to the given offset radius is 2.12% and 2.49% for δ = 10% and 5% respectively.

Dual Contouring has no guarantees for reconstructing the correct topology of a given implicit surface. This is a fundamental issue of all reconstruction methods based on structured volumetric discretization and is subject to many studies but has not been resolved entirely yet [ZGG22; RSA16; Gro16]. In our case, this means that the topology might be incorrect at concave creases with an angle that is too small to be resolved by the octree. As remeshing expects that the given offset mesh has the correct topology, it may create self-intersections in such areas. For δ = 10%, self-intersections appear in 2.9%, and for δ = 5% in 6.9% of all meshes. The discretization error remains small even in those regions because self-intersecting regions are pushed towards the concave creases. The appearance of self-intersections can be effectively reduced using the CGAL function remove_self_intersections from the CGAL component Polygon Mesh Processing [START_REF] Loriot | Polygon Mesh Processing[END_REF] to 0.19% and 0.65% for δ = 10% and 5% respectively, without impact on the first three significant digits of σm and εm. An example is depicted in Figure 8.

We also measure runtimes on Thingi10k, see Figure 9. For δ = 10%, we observe a correlation between input complexity and runtime for the topology-adapted Dual Contouring. In contrast, remeshing does not show any correlation, not even to the output complexity. The same behavior is observed for δ = 5%. The mean full runtime of our method is 91 and 148 seconds for δ = 10% and 5% respectively. While these values are insufficient to show correlation, they show that runtime tends to increase for smaller offset radii. The increased runtime for δ = 5% can be explained by the increasing offset complexity that needs to be captured. A larger offset radius is more likely to seal features like tunnels. The mean runtime is dominated by a few offsets that are compute intensive. Over 70% of all offset meshes are computed in less than 1 and 2 minutes respectively.

Analysis

In this section, we give an in-depth analysis of the two main parts, Dual Contouring on the topology-adapted octree, Section 6.2.1, and offset remeshing, Section 6.2.2.

Topology-Adapted Dual Contouring

The major added value of our octree in comparison to those that always refine to the maximum level, e.g. [START_REF] Pavi Ć | High-resolution volumetric computation of offset surfaces with feature preservation[END_REF], is that we avoid over-refinement in topologically simple regions. This is especially important when the offset has small features like narrow tunnels. This fact is exemplified by the couplingdown model that contains several drill holes. We generate the offset surface twice, once with the regular octree as described in Section 4 and once with the octree where the maximum level of subdivision is enforced by ignoring the disk criterion defined in Section 4.2. For a relative offset of 3%, an octree depth of 8 is required to capture all the tunnels within the drill holes, see Figure 10. In both scenarios, the topology was reconstructed correctly. However, generating the Dual Contouring mesh requires 34 seconds without our disk criterion but only 11 seconds with it. Additionally, the Dual Contouring mesh is substantially more complex and contains 523, 784 triangles, while the disk criterion reduces the complexity to 45, 392 triangles. After remeshing, which also takes longer without the criterion, 155 seconds vs. 110 seconds, both meshes have a low εm of 0.11% and 0.18%. The meshes then contain 104, 724 and 67, 840 triangles. Thus, remeshing can reduce the over-refinement caused by the octree very well but the mesh generated using the disk criterion remains less complex.

The advantage of the adaptive octree depends on the complexity of the offset topology and may not be that severe all the time. However, by using the adaptive octree, the user does not require a submitted to Eurographics Symposium on Geometry Processing (2023) priori knowledge of the offset topology and does not need to run the method multiple times until the octree depth is large enough to cover all topological features. Thus, this not only improves performance but also simplifies the method for the user.

Remeshing

The remeshing step is designed to reduce the distance and normal deviation from the offset meshes to the exact offset. Figure 11 illustrates remeshing at work on the anchor mesh with δ = 2% and σmax = 3 • . The discretization error εm is effectively reduced and converges after only 6 iterations.

We perform offsetting on the anchor model with offset radii δ = 15%, 10%, and 5%. For each radius, we set the maximum normal deviation to σmax = 20 • , 10 • , and 5 • . Both, σm and εm correlate to σmax, independent of the offset radius. The output meshes are depicted in Figure 12.

Applying Remeshing to Other Offset Methods

Mesh sharpening. Remeshing is completely independent of the octree constructed in Section 4 and can be applied to other methods that generate offset meshes. Alpha wrapping [PRH*22] is such a method: it aims to construct a manifold englobing volume around an input and for this generates an offset mesh with a user-defined precision. Although the method guarantees the absence of intersections with the input mesh, it does not recover sharp features in the output offset surface. Our remeshing step is relevant for postprocessing such an output. To avoid over-refinement in convex regions, we set a maximum normal deviation of 15 degrees during remeshing. For the degree-three corner model (Figure 13a) σm reduces from 6.5 • to 2.6 • and εm from 14% to 0.7%. A similar result is achieved on the screw nut, Figure 13d, with σm reducing from 10 • to 2.8 • and εm from 13% to 0.9%. In this process, the concave creases in the offset are reconstructed, as depicted in Figures 13c and13f.

Minkowski sums. Offsets can be also computed via Minkowski sums. Such a Minkowski-based approach was used by Campen and Kobbelt [START_REF] Campen | Polygonal boundary evaluation of Minkowski sums and swept volumes[END_REF]. We compare with this approach and show the effectiveness of our remeshing step by applying it to its output meshes, see Figure 14. On the filigree model, the Minkowski sum approximation consists of 732, 230 triangles for an absolute offset of approximately 0.07. For an absolute offset of 0.035 it consists of 1,051,354 triangles (Figure 14a). In comparison, our method yields 20,090 and 65,528 triangles, respectively (Figure 14c). When applying our remeshing approach to the Minkowski sum approximations, the complexity reduces to 22,892 and 73,772 triangles (Figure 14b). Additionally, triangle quality is significantly improved. While the Minkowski sum approximations contain selfintersections, ours contain none. Our remeshing successfully removes them also from the Minkowski sum approximations.

Limitations

The proposed approach guarantees that the output meshes are closed and combinatorially 2-manifold, but cannot guarantee the absence of self-intersections that are caused by offset meshes that are topologically different from the exact offset. This is an issue that all methods based on Dual Contouring or Marching Cubes share. In most cases, self-intersections could be resolved with a simple post processing. Self-intersections can be an issue for downstream applications and need to be resolved in the future. The user parameters provide a means to control the complexity of the output meshes, but a strict upper bound on the Hausdorff error is still an open problem. The remeshing step is effective at reconstructing sharp features, but its running time is hard to predict as the number and size of features present in the offset surface are unknown a priori.

Conclusion and Future Work

This paper introduced an approach for generating isotropic offset meshes from input triangle meshes. A first step applies the Dual Contouring method on an adaptive octree refined only where the offset topology is complex. A second step operates a novel remeshing method tailored to the offset that reconstructs sharp features and reduces the deviation between the normals of the output mesh and the normals of the offset while improving the shape of the mesh elements. The result is a reliable algorithm that operates even on defect-laden inputs, guarantees closed and combinatorially 2-manifold output meshes, and generates faithful offset meshes submitted to Eurographics Symposium on Geometry Processing (2023) with low approximation errors and well-shaped triangles. The novel offset-aware remeshing approach is also applied with success to offset meshes generated by other approaches.

In future work, we plan to explore a variant that generates anisotropic meshes with improved complexity-distortion tradeoffs. We also wish to explore alternatives to the axis-aligned octree data structure, such as unstructured tetrahedral meshes or binary space partitions.

  , the reliable generation of such a diagram in 3D is still open and only partial or approximate methods exist. In 3D, Boada et al. approximate the generalized Voronoi diagram using a voxel grid [BCMA08]. Hemmer et al. compute the exact 3D Voronoi diagram but only for lines [HSH10]. The complexity of computing the generalized 3D Voronoi diagram in Euclidean space is well described by Yap et al. [YSL12]. To our knowledge, there is no method that can compute either the exact diagram or an approximation with the precision required by our problem. Exact local offsets. Instead of using the complete generalized 3D Voronoi diagram, Aubry et al. use its local subset, namely the generalized spherical Voronoi diagrams around vertices, to compute the offset [ADMK17]

Figure 3 :

 3 Figure3: Our method handles triangles soups well and guarantees that the output is closed and manifold, even in complex scenarios.

Figure 4 :

 4 Figure 4: Closed meshes like caesar can have positive or negative offsets.

Figure 5 :

 5 Figure 5: We project the center of the cell (red) on each input triangle (black). If the spheres with offset radius centered at the projection points have a non-empty intersection (light red), the union of triangle offsets is guaranteed to be a star domain (blue).

  submitted to Eurographics Symposium on Geometry Processing (2023) D. Zint & N. Maruani & P. Alliez & M. Rouxel-Labbé / Feature-Preserving Offset Mesh Generation from Topology-Adapted Octrees

Figure 6 :

 6 Figure 6: Showcase of our results on Thingi10k.

Figure 7 :

 7 Figure 7: εm and σm for the models of Thingi10k sorted by εm.

Figure 9 :

 9 Figure 9: Runtime of our method measured on Thingi10k for δ = 10%.

  (a) Input (purple) and offset (blue) (b) Octree on maximum depth runtime: 190 seconds F = 105k εm = 0.11% σm = 1.07 • (c) Topology-adapted octree runtime: 121 seconds F = 68k εm = 0.18% σm = 1.31 •

Figure 10 :

 10 Figure 10: A hole in the input creates an arbitrarily thin tunnel in the offset surface. Our octree is automatically refined to fit small details on the offset.

Figure 11 :

 11 Figure 11: Remeshing reduces both, σm and εm. Renderings visualize εm after 0, 1, 2, and 9 remeshing steps.

Figure 12 :

 12 Figure 12: Anchor: δ = 15%, 10%, 5% (top to bottom), σmax = 20 • , 10 • , 5 • (left to right)

  εm = 0.23% σm = 2.3 •

Figure 14 :

 14 Figure 14: Comparison of our method to the approximated Minkowski sum (MS) from [CK10b]. Applying remeshing greatly reduces mesh complexity while keeping εm and σm low. It also produces better shaped triangles with more consistent angles.

submitted to Eurographics Symposium on Geometry Processing (2023)

Acknowledgments

This work has been co-funded by GeometryFactory and the national France Relance recovery plan. The work of Pierre Alliez is supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.