S Iseni

A Schmidt-Bleker

J Winter

K-D Weltmann

S Reuter

Atmospheric pressure streamer follows the turbulent argon air boundary in a MHz argon plasma jet investigated by OH-tracer PLIF spectroscopy

Keywords: atmospheric pressure plasma jet, planar laser-induced fluorescence, discharge propagation, plasma medicine, plasma jet effluent regime

An open question in the research of the dynamics of non-equilibrium cold atmospheric pressure plasma jets is the influence of ambient species on streamer propagation pathways. In the present work, by means of planar laser-induced fluorescence (PLIF), an atmospheric pressure argon plasma jet is investigated in a laminar and turbulent gas flow regime. The flow pattern is investigated with plasma on and plasma off. It is shown that in turbulent mode, the streamer path changes according to the flow pattern changes and the resulting changes in air abundance. From a comparison of an analytical diffusion calculation and LIF measurements, the air impurity boundary is determined. Most importantly, the origin of the streamer pathway is investigated in detail, by recording the flow pattern from OH-PLIF measurements and simultaneously measuring the streamer path by an overlay technique through emission measurements. It is shown that the streamer path is correlated to the turbulent flow pattern.

Introduction

Only a few years ago, so-called plasma bullet [START_REF] Laroussi | The plasma pencil: a source of hypersonic cold plasma bullets for biomedical applications[END_REF][START_REF] Teschke | High-speed photographs of a dielectric barrier atmospheric pressure plasma jet[END_REF] and guidedstreamer [START_REF] Robert | Characterization of pulsed atmospheric-pressure plasma streams (PAPS) generated by a plasma gun Plasma Sources[END_REF] jets came to be known. Recently, researchers have identified the fundamental processes behind these dynamics (see e.g. [START_REF] Xiong | Atmospheric pressure ionization waves propagating through a flexible high aspect ratio capillary channel and impinging upon a target Plasma Sources[END_REF][START_REF] Niemeyer | Wiley Encyclopedia of Electrical and Electronics Engineering[END_REF]). However, their propagation in air has continued to raise questions. In this work, the connection of gas flux regime, ambient species and streamer head propagation in a MHz atmospheric pressure argon plasma jet is investigated by laser spectroscopic methods. The diagnostic technique is based on planar laser-induced fluorescence (PLIF) on hydroxide molecules (OH) generated by the plasma jet itself due to water impurities. With stereoscopic imaging, both flow structure and streamer propagation are investigated simultaneously for a single discharge event. The study reveals that the flux field with plasma differs strongly from the flux field without plasma. Furthermore, it is most prominently and for the first time shown that the streamer propagates in the noble gas channel of a turbulent plasma jet effluent itself and follows the boundaries of ambient species impurities. Especially revealing are the differences observed for the laminar and turbulent case.

Experimental setup

In this work, the flux and streamer propagation of a MHz atmospheric pressure plasma jet (kINPen, neoplas tools GmbH, Germany) [START_REF] Weltmann | Atmospheric pressure plasma jet for medical therapy: plasma parameters and risk estimation Contrib[END_REF][START_REF] Reuter | Detection of ozone in a MHz argon plasma bullet jet Plasma Sources[END_REF] have been investigated. The feed gas is argon (Ar purity 99.999%), which may contain water impurities of up to 20 ppm originating from the gas pipe as reported in [START_REF] Winter | Feed gas humidity: a vital parameter affecting a cold atmospheric-pressure plasma jet and plasma-treated human skin cells[END_REF]. The OH PLIF measurements are performed with a frequency doubled dye laser (Cobra Stretch, Sirah Laser und Plasmatechnik GmbH, Germany) pumped with a 10 Hz pulsed Nd : Yag laser (Spectra-Physics, Inc.). The pulse duration is measured to be 8.1 ns after frequency doubling at 283 nm. The laser linewidth has been recently measured on a similar dye laser from the same company to be 1.4 pm [START_REF] Verreycken | Absolute calibration of OH density in a nanosecond pulsed plasma filament in atmospheric pressure He-H 2 O: comparison of independent calibration methods[END_REF]. Although this parameter remains crucial for the determination of absolute density [START_REF] Verreycken | Absolute OH density measurements in the effluent of a cold atmospheric-pressure Ar-H 2 O RF plasma jet in air Plasma Sources[END_REF], it does not influence qualitative measurements. The P1(4) transition of the OH(A, X)(1,0) at 283.5 nm is chosen as the energy transition due to its good absorption efficiency. Additionally, it is well separated from other lines. The linearity of the OH fluorescence signal versus the laser energy has been checked. The laser beam is expanded to a sheet via a cylindrical quartz lens and guided through the plasma jet effluent. The laser energy integrated along the 12 mm high sheet is measured to be 13 mJ. An aperture is mounted in front of the jet in order to illuminate only the area of interest. A CCD camera Imager Pro X (LaVision GmbH, Germany), combined with an external gated image intensifier and a 100 mm f/2.8 UV lens, is used as detection system. The fluorescence emission is observed at a 90 • angle to the p-polarized laser beam in order to reduce unwanted signals due to Rayleigh scattering. Additionally, a stereoscopic optical imaging system is coupled to the lens in order to record two images with each laser shot. By installing different filters for the two images respectively, different information can be gathered simultaneously. Here, on the first optical input, a BK7 neutral density filter is mounted, which allows observation from 350 to 1100 nm, where most of the plasma emission is found due to argon lines [START_REF] Winter | Feed gas humidity: a vital parameter affecting a cold atmospheric-pressure plasma jet and plasma-treated human skin cells[END_REF][START_REF] Reuter | From RONS to ROS: tailoring plasma jet treatment of skin cells[END_REF]. On the second optical input, a 10 nm FWHM band-pass filter centered at 313 nm is used to collect the fluorescence only of OH. Both signals are recorded on a vertical split CCD chip, so that the exposure time is identical for both signals.

The plasma jet frequency is phase locked with the laser pulse as well as the detection system, which allows phase-resolved optical emission spectroscopy (PR-OES). The exposure time is set to 200 ns, which is sufficient to obtain enough intensity for single shot acquisition. The optical aberrations such as distortion produced by the stereoscopic imaging system are corrected by Davis software [12] after the calibration procedure described in the manual [13]. The inhomogeneity of the laser sheet is measured by Rayleigh scattering and corrected by intensity normalization.

The observation of the flow pattern with plasma off is performed by PLIF using acetone as a tracer. As described in [START_REF] Lozano | Acetone: a tracer for concentration measurements in gaseous flows by planar laser-induced fluorescence[END_REF], acetone fluorescence is a suitable tracer to image turbulent flows. During the measurements, the feed gas is bubbled through acetone into the plasma device. The acetone molecules are excited by the laser beam at 283.5 nm, and the fluorescence from 350 to 800 nm is collected by an intensified CCD camera and recorded. Observing the flow pattern with plasma-on requires the use of a different tracer, since acetone would be destroyed by the reactive species in the plasma. In the past, hydroxyl fluorescence was already used in combustion with good success on turbulence studies [START_REF] Tanahashi | Simultaneous CH-OH PLIF and stereoscopic PIV measurements of turbulent premixed flames[END_REF]. Moreover, one report stated that the hydroxyl molecule lifetime is 1 ms to a few milliseconds at atmospheric conditions [START_REF] Kanazawa | LIF imaging of OH radicals in DC positive streamer coronas[END_REF][START_REF] Pei | A study on the temporally and spatially resolved OH radical distribution of a room-temperature atmospheric-pressure plasma jet by laser-induced fluorescence imaging[END_REF]. The hydroxyl radical is produced by the plasma itself from water impurities within the feed gas [START_REF] Winter | Feed gas humidity: a vital parameter affecting a cold atmospheric-pressure plasma jet and plasma-treated human skin cells[END_REF], meaning that the plasma generates the tracer itself and no external species need to be admixed for the flow measurements. Recently, hydroxyl molecule densities on-axis were measured by UV absorption on the same discharge [START_REF] Verreycken | Absolute OH density measurements in the effluent of a cold atmospheric-pressure Ar-H2O RF plasma jet in air Plasma Sources[END_REF]. The results show that the density remains constant up to a 7 mm distance from the nozzle. It can thus be assumed that the hydroxyl density does not decay significantly along the first 10 mm from the nozzle [START_REF] Verreycken | Absolute OH density measurements in the effluent of a cold atmospheric-pressure Ar-H2O RF plasma jet in air Plasma Sources[END_REF]. For these reasons, the hydroxyl radical, for the present studies, can be used as a fluorescent tracer.

LIF on acetone-plasma off measurements

Figure 1 shows the fluorescence of acetone for three different feed gas flow rates. As can be seen, at both 0.5 and 1 standard liter per minute (slm) the flow is laminar, while at 3 slm the flow is turbulent. The indicated Reynolds numbers Re = 2r 0 v av /ν are defined via the capillary radius r 0 and the kinematic viscosity ν. The critical Reynolds number (Re) for the transition from the laminar to turbulent regime in pipe flows is Re C = 2300, which is consistent with the observed flow regimes.

OH fluorescence and quenching by ambient air-plasma on measurements

The radiative lifetime of the OH(A, X)(1,0) is known to be around τ = A -1 = 748 ns, where A is the Einstein coefficient for spontaneous emission. The emission is strongly quenched by atmospheric species [START_REF] Verreycken | Absolute OH density measurements in the effluent of a cold atmospheric-pressure Ar-H2O RF plasma jet in air Plasma Sources[END_REF]: the combined quenching rate (Q) and vibrational energy transfer (VET) coefficient of OH(A, X)(1,0) by argon (0.27 × 10 -11 cm 3 s -1 [START_REF] Dilecce | Laser diagnostics of high-pressure discharges: laser induced fluorescence detection of OH in He/Ar-H 2 O dielectric barrier discharges Plasma[END_REF]) is two orders of magnitude lower than for nitrogen (k [START_REF] Williams | Collisional vibrational energy transfer of OH (A 2 + , v = 1)[END_REF]). Despite that, the quenching rate of OH fluorescence by water is in the same order of magnitude as dry air, and ambient humidity has a negligible effect on the OH fluorescence quenching due to a much lower density compared to air. In consequence, OH(A, X)(1,0) has a longer lifetime when it is mixed with argon than in presence of air. The OH originates dominantly from the dissociation of water impurities in the feed gas.

N 2 Q = 23.3 × 10 -11 cm 3 s -1 [20]), oxygen (k O 2 Q = 20.6×10 -11 cm 3 s -1 [20]) and water (k H 2 O Q = 66 × 10 -11 cm 3 s -1
Even at what is considered to be dry feed gas conditions, the quantity of water molecules lies in the order of some 10 ppm. When the gas leaves the jet, the OH concentration is distributed relatively homogeneously along the diameter of the nozzle (see figure 2(c)). Similar observations have been reported earlier. Here the main hydroxyl production is localized near the powered electrode in case of water impurities in the feed gas [START_REF] Yonemori | Measurement of OH density and air-helium mixture ratio in an atmospheric-pressure helium plasma jet[END_REF]. Additionally, the hydroxyl density remains presumably constant until several mm distance [START_REF] Verreycken | Absolute OH density measurements in the effluent of a cold atmospheric-pressure Ar-H2O RF plasma jet in air Plasma Sources[END_REF]. The spatial distribution can be expected to resemble the acetone distribution as shown above (see figure 1). From these as well as the quencher considerations, it can be safely assumed that the LIF signal reveals rather the influence of the quenching by the ambient species than the ground state OH molecule distribution. The measured data therefore reflects the argon flow pattern and the argon air boundary shape, as will be described in the following.

The diffusion of ambient species into the effluent of the plasma jet can be estimated by solving the convection diffusion equation v(r)∂ z n A + D r n A = 0 with the boundary condition n A = θ(r -r0)n 0 A for the ambient species density, where θ(x) is the Heaviside step function θ(x) = 0 for x < 0, θ(x) = 1 for x 0. Here, D is the diffusion coefficient and v(r) a velocity field with a radial dependence and r the radial part of the Laplace operator. The solution in the case of a constant velocity profile is given in [START_REF] Saìnchez | Exact solutions for transient mixing of two gases of different densities[END_REF] and can be evaluated using an infinite series. For a laminar jet flow, the velocity field in the near field of the jet can be approximated by a Poiseuille flow profilev(r) = 2v av (1r 2 /r 2 0 ). In order to account for this flow profile, non-dispersive path mapping [START_REF] Schmidt-Bleker | Non-dispersive path mapping approximation for the analysis of ambient species diffusion in the near-field, near-axis region of laminar jets in preparation[END_REF] can be applied, mapping the constant velocity solution according to the flow field. The contour lines of the resulting ambient air density n A /n 0 A for an average velocity of v av = 4.1 m s -1 and diffusion coefficient D = 0.2 cm 2 s -1 are shown in figure 2(a). Assuming a constant OH-density in the jet's effluent, the fluorescence intensity of OH(A, X)(1,0) can be estimated to be proportional to

I LIF ∼A/(A + n O 2 k O 2 Q + n N 2 k N 2 Q + n H 2 O k H 2 O Q )
as depicted in figure 2(b). In figures 2(c)-(e), the single shot LIF signal at flow rates of 0.5, 1 and 1.5 slm are displayed. The signal is a pure LIF signal, as spontaneous emission from excited hydroxyl OH(A, X)(1,0) is negligible. The analytically obtained triangular emission pattern (figure 2(b)) agrees very well with the experimentally observed OH emission at a gas flux of 0.5 slm (figure 2(c)). A sharp edge in the emission pattern occurs at an ambient air density of approximately 0.01-1%. This shows that due to the strong hydroxyl fluorescence quenching by air, the contour of the argon-air boundary can be determined and reflects the relative space-resolved argon density. The plasma does not change the flow regime at this flow rate. For a flow rate of 1 slm, the occurrence of Kelvin-Helmholtz (KH) instabilities is observed, leading to a breaking of the laminar jet downstream. KH instabilities are likely to appear at the mixing layer of two fluids with different densities and/or in a shear flow [START_REF] Kull | Theory of the Rayleigh-Taylor instability[END_REF]. Both mechanisms may be triggered when heating of the gas by the plasma occurs: its density is decreased, while its expansion leads to an increased shear.

Electro hydrodynamic forces can also trigger flow regime transitions. In [START_REF] Sarron | Plasma plume length characterization 21st Int[END_REF] it is observed that the flow regime of a helium jet is significantly affected by the plasma while the temperature remains at around room temperature. It is noted that helium jets may be affected more strongly by equal electro hydrodynamic forces, as their momentum is ten times less compared to argon jets at identical flow conditions. It was recently observed that the air mixing in an argon plasma jet similar to the one used in the present study is significantly increased when the plasma is turned on, which agrees with our findings [START_REF] Van Gessel | Electron properties and air mixing in radio frequency driven argon plasma jets at atmospheric pressure[END_REF]. Their gas temperature is around 480 K. They also attribute this effect to the occurrence of KH instabilities and/or electro hydrodynamic forces.

At 3 slm, the flux seems to remain fully turbulent as in the plasma-off case.

Turbulent flow and streamer propagation

The observation of the turbulent behavior of the effluent flow is possible only if the acquisition time is much faster than the turbulence time scale. At 200 ns exposure time, the flow dynamics is static, but the discharge propagation is fast enough to grow and develop itself. Therefore looking at the effluent turbulence pattern as well as the streamer path within the same single time shot leads to fundamental information about the discharge propagation. Figure 3 shows a superimposed image of the streamer and the hydroxyl fluorescence. The hydroxyl fluorescence resembles the argon flow pattern. Both are following the same path and-especially at greater distances to the nozzle, where the argon abundance has the same dimension as the streamer thickness-the patterns are identical. On figure 3(b) and (c), which corresponds to 1 slm and 3 slm flow rates respectively, the flow patterns are not laminar and result in density gradients of argon which form a nonrectilinear argon channel. Also the streamer propagation does not follow a straight path. For the 3 slm case, it is clear that the turbulent flow has an impact on the streamer morphology (see figure 1(c)). This argon channel is about 2 mm in size at the nozzle and thins down as it propagates along the jet-axis. The experimental results show that the streamer propagates into the argon channel from the wider side and ends its propagation within the very thin argon channel surrounded by air impurities from diffusion, at the space where the argon density remains the highest. This means that the streamer is always bound by argon abundance, as shown in figure 3. In the case of 0.5 slm (figure 3(a)), the argon flow is laminar, inducing a rectilinear argon channel where the streamer propagates linearly. The argon flow is only clearly visible along the first 5 mm because of the induced fluorescence, which is quenched by the diffusion of ambient species as described previously. Moreover, we observe that the streamer position is well aligned on the axis of a right circular cone resulting from the gradient of the air diffusion into the argon channel. It is especially true towards greater distances from the nozzle. These observations are in perfect agreement with the fact that the streamer propagation follows the highest argon density which obviously is located on the axis of the channel. However, those measurements have been made over a 200 ns gate which means that the flow dynamics remains static with respect to the streamer head propagation. In order to confirm the results observed from single shot acquisitions, successive acquisitions of streamer propagations are recorded with the same gating time with an interval of 100 ms between two shots. At this acquisition frequency, the flow dynamics has time to develop from one shot to the other. Thus, it is ensured that between two consecutive acquisitions, the flow pattern must differ from one to the other in the case of turbulence or non-stationary flow.

Accumulation of 50 streamers is recorded for 0.5 slm and 3 slm flow rates, which correspond to a laminar and a turbulent flow, respectively. The results are presented in figure 4. In the turbulent case at 3 slm, the streamers' spatial distribution is spread over a wider area compared to the laminar case at 0.5 slm. This observation is even more pronounced after a 6 mm distance from the nozzle, where the structure of the turbulence is developing, as shown in figure 2(c). For the laminar case, (figure 4(a)), the position of the streamers remains constant. This is in agreement with the laminar flow pattern measured for the same conditions and presented in figure 2(a) . It thus becomes evident that the atmospheric streamers are bound by argon and propagate through this argon channel. A previous study on streamer propagation in argon flow surrounded by air has been published [START_REF] Hofmann | Transitions between and control of guided and branching streamers in dc nanosecond pulsed excited plasma jets[END_REF]. The authors report that the reason for the propagation of a guided streamer in argon in the ambient air is due to the interface between these two gases. They observe the brightest emission on the edge on the plasma effluent and identify it as an emission from the second positive system of nitrogen, N 2 (C), mainly resulting from the efficient energy transfer of the near-resonant transition between the argon metastables of the ( 3 P 0,2 ) states and N 2 (X) to produce N 2 (C) [START_REF] Nguyen | Rotational and vibrational distributions of N 2 (C 3 u) excited by state-selected Ar(3P2) and Ar(3P0) metastable atoms[END_REF]. In this work, N2(C) emission position differs from the referred study and is located in the center of the effluent. Moreover, the highest emission intensity was detected in the range from 4 mm up to 7 mm from the nozzle. This is in excellent agreement with figure 2(c) ), where one sees that OH quenching by ambient molecules reaches the on-axis. Therefore it can be assumed that metastables of argon are, in the present case, located within the argon channel where the streamer develops. One reliable reason for the difference observed to [START_REF] Hofmann | Transitions between and control of guided and branching streamers in dc nanosecond pulsed excited plasma jets[END_REF] is attributed to the presence of a grounded ring electrode in the plasma source which must induce a different electric field and thus different ionization mechanisms. Due to the high frequency in the present case, streamer events will benefit from previous ionization and excitation phenomena. The memory effect of metastables and ions contributes largely to the development of the discharge within the argon dominating channel. This is also valid for the turbulent cases as the vortexes remain constant on the time scale of streamers.

A simulation of discharge development in the case of a pure helium laminar flow surrounded by air without diffusion explains that socalled streamer bullet propagation is more favourable within the helium tube than in air because of the ionization coefficient differences of each gas [START_REF] Boeuf | Dynamics of a guided streamer ('plasma bullet') in a helium jet in air at atmospheric pressure[END_REF].

The similarity in streamer path and gas flow path was also found by [START_REF] Robert | Rare gas flow structuration in plasma jet experiments Plasma Sources[END_REF] in pulsed atmospheric pressure plasma streams operated in helium and neon. Interestingly, they found clear evidence that under their experimental conditions, the plasma determines the gas channel behavior and not, as commonly thought, vice versa [START_REF] Sarron | Plasma plume length characterization 21st Int[END_REF][START_REF] Robert | Rare gas flow structuration in plasma jet experiments Plasma Sources[END_REF][START_REF] Pouvesle | Plasma jets for biomedical applications: PAPS control of the rare gas jets 66th Gaseous Electronics Conf[END_REF]. Since the mass density of argon is ten times higher than for helium, this plasma-induced gas channel forming is not directly transferable to the argon plasma jet investigated in this work. Hence, for the laminarto-turbulent transition in the 1.0 slm case when the plasma is switched on, it is not clear whether the streamers are responsible for inducing the turbulence of the argon flow or whether the argon flow becomes turbulent itself by pure fluid dynamics. Experiments with the focus on this effect need to be performed in future. The studies presented in the present work clearly show that the path of highest argon density, lowest air impurity density and streamer propagation are identical.

Conclusion

In this study, it was shown by means of PLIF on an atmospheric pressure argon plasma jet that in turbulent mode, the streamer path changes according to the flow pattern changes and the resulting changes in air abundance. Furthermore, with known OH quenching coefficients from a comparison of analytical diffusion calculation and LIF measurements, the air impurity boundary can be determined. Comparing the plasma-on and plasma-off case, it was shown that in the laminar case, the development of KH instabilities is promoted, while in the turbulent flow, the flow patterns remain comparable. Most importantly, it was shown that the streamer pathway is defined by the ambient air diffusing into the effluent. From a convection diffusion model, a steep increase of ambient air particle density in radial direction was revealed. This gradient forms a boundary which leads to the observed streamer enclosure within the argon channel.

Figure 1 .

 1 Figure 1. Space resolved single shot fluorescence signal of acetone for different argon gas flow rates of (a) 0.5 slm/Re= 465, (b)1.0 slm/Re= 930 and (c) 3.0 slm/Re= 2790. The jet nozzle is located at an axial distance of zero. The plasma was switched off.

Figure 2 .

 2 Figure 2. Analytical estimation of the ambient species density (a) and the resulting analytical OH emission pattern (b). Space resolved fluorescence signal of hydroxyl for different argon gas flow rates of (c) 0.5 slm, (d) 1.0 slm and (e) 3.0 slm. The jet nozzle is located at an axial distance of zero. The plasma was switched on. Single shot acquisition.

Figure 3 .

 3 Figure 3. Overlay of the spatial OH fluorescence signal and of the spatial plasma emission (visible and near-infrared) for different argon gas flow rates of (a) 0.5 slm, (b) 1.0 slm and (c) 3.0 slm. Both signals were obtained simultaneously on the same detector.

Figure 4 .

 4 Figure 4. Spatial streamer distribution in the jet effluent for an argon gas flow rate of (a) 0.5 slm and (b) 3.0 slm. The number of accumulations was 50. (The slight vertical tilt is due to the alignment of the plasma jet.)

Acknowledgments

Discussions with Dr E Robert and Professor J-M Pouvesle in the frame of collaborative project between Leibniz Institute for Plasma Science and Technology INP Greifswald e.V. and Groupe de Recherches sur l'Energétique des Milieux Ionisés (GREMI) are gratefully acknowledged. The authors gratefully acknowledge the funding by the Federal German Ministry of Education and Research (BMBF, grant number 03Z2DN12).

This is the Accepted Manuscript version -post-print-of an article accepted for publication in Journal of Physics D: Applied Physics. This article may be downloaded for personal use only. Any other use requires prior permission of the authors and IOP Publishing. This article may be found at doi:10.1088/0022-3727/47/15/152001. IOP Publishing is not responsible for any errors or omissions in this version of the manuscript or any version derived from it.


David Desobry

Yoann Coudert-Osmont

Etienne Corman

Nicolas Ray

Dmitry Sokolov

Designing 2D and 3D Non-Orthogonal Frame Fields

Keywords: Geometry processing, vector field design, quadrilateral remeshing, hexahedral remeshing

We present a method for direction field design on surface and volumetric meshes supporting non-orthogonality. Our approach is a generalization of the representation of 3D cross fields in spherical harmonic basis. As such it induces a geometrically meaningful measure of smoothness, allows orthogonality control by a simple parameter and enables orientation constraints of a single direction. To the best of our knowledge this is the first work to propose non-orthogonal 3D frame field design. We demonstrate the applicability of our method to generate anisotropic quadrangular and hexahedral meshes which are particularly useful for remeshing CAD models.

Introduction

Designing tangent vector fields is a fundamental component of geometry processing. If it is predominantly used for surface remeshing [START_REF] Kälberer | Quadcover-surface parameterization using branched coverings[END_REF] it has also applications in stylized rendering [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF], BRDF modification [START_REF] Brandt | Modeling nsymmetry vector fields using higher-order energies[END_REF], fabrication [START_REF] Sageman-Furnas | Chebyshev nets from commuting polyvector fields[END_REF] and sketch based modeling [START_REF] Iarussi | Bendfields: Regularized curvature fields from rough concept sketches[END_REF][START_REF] Bessmeltsev | Vectorization of line drawings via polyvector fields[END_REF]. Vector field design is often a single step of a more involved pipeline. Since each application has its own needs it is important to provide robust and flexible algorithms.

Typical applications like parametrization require to design N-Directions vector fields: 2D frame fields with 2-directions for quad-meshes [START_REF] Kälberer | Quadcover-surface parameterization using branched coverings[END_REF] and 3D frame fields with 3-directions on a volumetric domain for hex-meshing [START_REF] Nieser | Cubecover-parameterization of 3d volumes[END_REF]. Early works only allowed the design of orthogonal fields, often called cross fields, to guide the orientation of the quads or hexes. However, the orthogonality constraint is often too restrictive and does not offer enough design freedom.

Meshing of CAD models is a typical example of how limiting orthogonality is. In these models it is common to encounter chamfers and sharp corners. Cross fields simply cannot meet these design constraints as demonstrated in Figure 4 and Figure [START_REF] Ray | N-symmetry direction field design[END_REF]. Many articles propose custom solutions to stray away from the dull isotropic-right-angled frames, by creating a new representation [START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF], by deforming the mesh [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF] or by modifying the metric tensor [START_REF] Jiang | Frame field generation through metric customization[END_REF]. When doing so, three challenges must be overcome.

First they must design a new representation or framework in which to define their vector fields. This representation is highly impacted by the dimensionality of the targeted space. In many cases, representation well-suited for surfaces, like the representation vector for 2D cross field with complex numbers [START_REF] Ray | N-symmetry direction field design[END_REF][START_REF] Knöppel | Globally optimal direction fields[END_REF][START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF], cannot be directly extended to volumetric meshes. Instead a completely different approach was introduced taking advantages of the symmetries of the spherical harmonic basis [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF]. This situation leads to large body of literature declining vector field design in various setups: surface or volumetric meshes, isotropic or anisotropic, orthogonal or non-orthogonal, sometimes cube symmetry or higher order symmetry. . . with no or little overlaps in these methods.

A common design requirement is to find the smoothest vector field interpolating between constraints. So once the representation is fixed the second challenge is to build a geometrically meaningful measure of smoothness of the field. This is of key importance as a non-geometrically based measure might lead to non-intuitive vector field design.

The third challenge is to be able to impose user defined constraints within the chosen representation. Such constraints may be as-simple-as a hard constraint on the vector field, more subtle like fixing one direction of the frame and let the other directions rotate freely or allow control over the orthogonality of the frame.

Contributions. In this paper, we represent sets of directions (or direction fields) as a polynomial function on a circle or a sphere (depending on the dimensionality). This representation is directly inspired by the spherical harmonic representation of 3D cross fields introduced in [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF]. By designing polynomials, we can choose the number of symmetries and the shape of the frame as illustrated in Figures 2 and8. This representation allows rotation of each direction independently which creates a parametric space of anisotropic non-orthogonal frames in any dimension. Unlike previous methods, ours addresses efficiently the three challenges:

1. The field representation is readily working for surfaces and volumetric meshes; 2. It offers a geometrically meaningful notion of smoothness; 3. Any type of direction fields can be designed and a simple parameter allows control over the orthogonality.

The rest of the paper will be organized as follow. Section 1 gives a general "recipe" to design direction fields of any dimension, thanks to a representation based on a sum of atomic polynomials. This theoretical framework is then applied to 2D frame field design in Section 2 and to generation of 3D frame fields in Section 3. In these sections, we will demonstrate the usefulness of our representation for quad-meshing and hex-meshing.

Related work

Frame fields on surfaces and volumes have a wide range of applications in computer graphics. In this paper, we will focus on global parametrization. Cross fields and more generally frame fields have been largely used in two dimension for quadrangular remeshing [START_REF] Kälberer | Quadcover-surface parameterization using branched coverings[END_REF], design of 6-symmetry vector fields can produce highly regular triangular remeshes [START_REF] Nieser | Hexagonal global parameterization of arbitrary surfaces[END_REF] and 3D frame fields are a common tool for hexahedral mesh generation [START_REF] Nieser | Cubecover-parameterization of 3d volumes[END_REF].

Although symmetric vector fields of two and three dimensions share many applications, their representation are vastly divergent. In fact, many well established 2D vector field design methods cannot be extended to 3D.

Frame field design on surfaces

Cross fields have been explored extensively in the past decade. They were first introduced in computer graphics to place hatches in non photorealistic rendering [START_REF] Hertzmann | Illustrating smooth surfaces[END_REF]. They were soon generalized to n-symmetry vector fields [START_REF] Palacios | Rotational symmetry field design on surfaces[END_REF][START_REF] Ray | N-symmetry direction field design[END_REF] by representing a cross as a n th root of a complex number. The smoothness of a field is simply obtained by measuring the similarity of pairs of adjacent frames. This leads to a simple quadratic optimization problem allowing alignment constraints [START_REF] Ray | N-symmetry direction field design[END_REF][START_REF] Knöppel | Globally optimal direction fields[END_REF] or singularity constraints [START_REF] Crane | Trivial connections on discrete surfaces[END_REF]. Cross field design surveys can be found in [START_REF] Vaxman | Directional field synthesis, design, and processing[END_REF] and [START_REF] De Goes | Vector field processing on triangle meshes[END_REF].

Liu et al. [START_REF] Liu | General planar quadrilateral mesh design using conjugate direction field[END_REF] proposed a generalization of cross fields by allowing them to be non-orthogonal. The goal was to generate quad meshes with the possibility for each element to stray away from a perfect square. However, this method requires to compute an explicit matching between the vector set to disambiguate the symmetry. More recent works generate anisotropic and nonorthogonal frame fields by deforming the mesh or designing a new metric [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF][START_REF] Jiang | Frame field generation through metric customization[END_REF]. Processing metric tensors is a costly which implies to solve an intricate optimization problem with multiple objectives. These methods could be generalized to volumetric meshes but at great computational cost since symmetric 3-tensors have twice as much degrees of freedom as 2-tensors.

Closer to our work is the representation of n-symmetry, nonorthogonal vectors as a root of a complex polynomial [START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF][START_REF] Diamanti | Integrable polyvector fields[END_REF]. The authors remarked that a set of 2D directions can be seen as the set of roots of a complex polynomial. In this setup, a deformation of the frame amounts for a modification of the polynomial coefficients. Therefore, they use polynomial coefficients to uniquely represent each frame. Like the representation described in the present paper, the complex polynomial representation is invariant by reordering of the direction set. However, this choice of variables blurs the link with the geometry of the frame and create two limitations not present in our work. First, in order to find smooth frame fields, a smoothness energy is applied to the polynomial coefficients. Although the smoothness of the field implies the smoothness of the coefficients, this notion of smoothness is not geometrically rooted and may introduce bias. An in depth discussion is provided in Section 2.5. Second, this representation does not generalize to higher dimensions.

3D frame fields and volumetric representation

The 3D cross fields, sometimes referred to as octahedral fields, are commonly used for hexahedral meshing [START_REF] Nieser | Hexagonal global parameterization of arbitrary surfaces[END_REF]. A single frame consists of three orthogonal vectors and their opposite, thus a 3D cross has the same symmetries than a cube. Huang et al. [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF] proposed to represent each individual frame by the spherical function (x, y, z) → x 4 + y 4 + z 4 whose level sets are cube-like. While seemingly different than the two dimensional case, Ray et al. [START_REF] Ray | Practical 3d frame field generation[END_REF] noted that it can be understood has a generalization of 2D complex representation vector. In this paper, we will build upon this observation.

As for the 2D case, generalization to n-symmetry 3D vector fields [START_REF] Shen | Harmonic functions for rotational symmetry vector fields[END_REF] and to anisotropic orthogonal frame field [START_REF] Chemin | Representing three-dimensional cross fields using 4th order tensors[END_REF][START_REF] Palmer | Algebraic representations for volumetric frame fields[END_REF] have been proposed by studying the properties of the spherical harmonic basis. However, to the best of our knowledge, there is no existing work on non-orthogonal 3D frame field design. Moreover, a representation unifying 2D and 3D direction fields allowing independent direction constraints is still missing.

Direction field of any dimension as spherical polynomials

Given a simplicial n-complex F , we consider each facet f ∈ F as a discrete tangent space. Consequently, in this paper tangent vector fields will be piecewise constant. To be more specific, on a triangle mesh our tangent vector fields are defined on faces of triangles and on a tetrahedral mesh tangent vector fields are constant per tetrahedron. This is a standard setup for global parametrization and structured meshing [START_REF] Kälberer | Quadcover-surface parameterization using branched coverings[END_REF][START_REF] Nieser | Cubecover-parameterization of 3d volumes[END_REF]. We denote E the set of pairs of adjacent facets; namely i, j ∈ F are adjacent if and only if (i, j) ∈ E. For imposing boundary conditions, we will need ∂F , the set of boundary facets.

To ease the exposition, we restrict ourselves to flat domains of R n . All the concepts introduced in this paper can be transposed for curved surfaces. One only has to define a notion of parallel transport as done in previous papers [START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF][START_REF] Ray | N-symmetry direction field design[END_REF].

Representation of a single frame

We will call a N-Direction set an unordered set of N unit tangent vectors and their opposite; namely

{u 1 , -u 1 , u 2 , -u 2 , . . . , u N , -u N }.
A N-Direction field is an assignment of a direction set per facet. Our goal is to find the smoothest N-Direction field on a triangle or tetrahedral mesh satisfying user-defined direction constrains. As a direction set is unordered, measuring the similarity of two sets can be rather painful as it requires to correctly match each set. Therefore, we need a representation of a single direction set that is invariant to sign change and reordering. Inspired by the 2D and 3D cross fields representation introduced in [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF][START_REF] Ray | Practical 3d frame field generation[END_REF], we will identify a direction to a spherical polynomial whose shape and symmetry is identical.

As a direction is a set of two opposite vectors, it must be represented by an even degree polynomial in order to match the symmetry. Indeed, given a unit vector u, the spherical polynomial p(s) = ⟨u, s⟩ 2k for k ≥ 1 integer is independent of the sign of u.

To represent the N-Direction set {u i , -u i | i = 1 . . . N}, we simply sum over these atomic polynomials:

p(s) = N i=1 ⟨u i , s⟩ 2k . (1)
Figure 1 illustrates the polynomial representing a unique direction in polar and spherical coordinates. Most interesting to parametrization applications are the 2D and 3D frames illustrated in Figures 2 and8 for k = 2.

Dealing with polar or spherical polynomial may seem difficult but it is made simple by decomposing it into canonical orthonormal bases.

Tangent space of a Direction set

Considering a N-Direction set for all possible unit vectors spans a polynomial subspace denoted P 2k,N :

P 2k,N :=        N i=1 ⟨u i , s⟩ 2k | ∀u i ∈ S n        .
The function space P 2k,N contains only homogeneous polynomials of degree 2k. As such, it is a finite dimensional space and is spanned by a finite number of orthogonal functions Y ℓ . Later, we will make use of the Fourier basis for 2D direction fields and the Spherical Harmonic basis for 3D direction fields. Therefore,

p(s) = L ℓ=1 c ℓ (U)Y ℓ (s).
Crucially, the coefficients c ℓ (U), like the polynomial p(s), are independent to reordering and sign changes of the vectors in the set U. They only depend on the positions of the vectors as a whole. This representation greatly simplifies the computation of a distance between two N-Direction sets as it suffices to compute a distance between their polynomial representations. The representation polynomials have the same shapes as the sets (see Figures 2 and8 for an example of frame field) so this distance has a clear geometrical meaning.

Using the fact that the basis Y ℓ is orthogonal, the distance between a polynomial p 1 representing the set of directions U 1 and p 2 representing the set of directions U 2 is simply the norm of the difference of the coefficients c:

∥p 1 -p 2 ∥ 2 2 = S n (p 1 (s) -p 2 (s)) 2 ds = L ℓ=1 (c ℓ (U 1 ) -c ℓ (U 2 )) 2 .

Piecewise constant polynomial PolyVector

A N-Direction field assigns a set U i per facet i ∈ F . Our representation replaces the sets U i by the coefficients c ℓ (U i ), ℓ = 1, . . . , L. Our goal is to design smooth N-Direction fields. A N-Direction field is considered perfectly smooth across an edge if its polynomial, and thus its coefficients c, are equal. Following this intuition, we introduce a smoothness energy, akin to a discrete Dirichlet energy, by accumulating the distance between adjacent frames:

E s (U) = (i j)∈E L ℓ=1 c ℓ (U i ) -c ℓ (U j ) 2 . (2)
Designing smooth direction fields interpolating constraints amounts for minimizing the smoothness energy E s with direction constraints directly imposed on the vector set U. Making design problems much more accessible than before.

Degree of the representation polynomial

The degree of the polynomial is a parameter that influences the shape of the frame and the complexity of the optimization problem. There is a lower bound for the polynomial degree. In [START_REF] Shen | Harmonic functions for rotational symmetry vector fields[END_REF], the authors studied spherical polynomial for representing N-Direction sets evenly distributed in 3D. They remarked that a polynomial of degree at least 4 was needed to represent 3 orthogonal directions and a 6 degree polynomial was necessary for 6 directions forming the vertices of a regular dodecahedron. Moreover, in order for the N directions of a N-Direction set to rotate independently the polynomial must be of degree high enough to represent these variations. However, the complexity of the coefficients c increases with the polynomial degree thus the minimization of the smoothness energy in Eq. 2 becomes more challenging .

In the rest of the paper, we focus on the special case of frame fields of 2 and 3 dimensions. We will choose the lowest acceptable polynomial degree, in the present case degree 4 or k = 2 in Eq.( 1). Figures 2b and8b shows that polynomials of degree 4 are sufficient to catch deformations of a 2D or 3D frame.

2D frame fields

In this section, we focus on frame fields (or 2-Direction fields) on 2D domains. First introduced in [START_REF] Panozzo | Frame fields: Anisotropic and non-orthogonal cross fields[END_REF], they are extremely useful for quad-meshing. For 2D domain the representation polynomial p is defined on a circle and the set of admissible tangent frames P 2k,2 is decomposed on the Fourier basis. We will see that with this choice of basis, and when the two directions are orthogonal, our representation is exactly equal to the complex representation of cross fields used in many articles [START_REF] Palacios | Rotational symmetry field design on surfaces[END_REF][START_REF] Ray | N-symmetry direction field design[END_REF][START_REF] Knöppel | Globally optimal direction fields[END_REF][START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF].

One direction as a polynomial of degree 4

As noted in [START_REF] Ray | Practical 3d frame field generation[END_REF], 2D cross fields can be represented by the polar polynomial of degree four: p(s) = s 4 1 + s 4 2 . In fact, this is the polynomial of the lowest degree presenting the same symmetry as a cross. For this reason, we will use degree 4 polynomials to represent our frame field.

A 2D frame is the set of two directions {u, -u, v, -v} parametrized respectively by the angles α and β, such that u = (cos α, sin α) and v = (cos β, sin β). These notations are summarized by Figure 2. Following Eq. (1) our representing polynomial of degree four reads as the sum:

p(s) = ⟨u, s⟩ 4 + ⟨v, s⟩ 4 . (3)
Figure 2 gives an illustration of the shape of an orthogonal frame (Fig. 2a) and a non-orthogonal frame as represented by the polar polynomial in Eq. (3).

Fourier basis

The representing polynomial in Eq. ( 3) is defined on a unit circle. It seems natural to use the Fourier basis as an orthogonal basis to decompose this periodic signal. We will note c ℓ,1 (α, β) the coefficients corresponding to the basis functions cos(ℓt) and c ℓ,2 (α, β) the coefficients corresponding to the basis functions sin(ℓt) where ℓ ∈ N is the frequency index.

The sine/cosine basis appears naturally when expressing s = (cos t, sin t) in the polar coordinate system. The direction polynomial in u simply reads:

⟨u, s⟩ 4 = |cos t cos α + sin t sin α| 4 = 3 4 + 1 2 cos(2t) cos(2α) + 1 2 sin(2t) sin(2α) + 1 8 cos(4t) cos(4α) + 1 8 sin(4t) sin(4α)
Thus the coefficients of the decomposition of polynomial frame boils down to:

c 2,1 (α, β) = 4 cos(2α) + 4 cos(2β) c 2,2 (α, β) = 4 sin(2α) + 4 sin(2β) c 4,1 (α, β) = cos(4α) + cos(4β) c 4,2 (α, β) = sin(4α) + sin(4β) (4)
Our frame is represented by only four coefficients in frequencies ℓ = 2 and ℓ = 4. Crucially, the frame is orthogonal (i.e. β = α + π/2) if and only if the coefficients of second order are zero. In this case, the frame is solely represented by the vector (cos(4α), sin(4α)), thus we recover the classic complex representation of cross fields [START_REF] Palacios | Rotational symmetry field design on surfaces[END_REF][START_REF] Ray | N-symmetry direction field design[END_REF].

Designing 2D frame fields for quad meshing

In this section, we will set up an optimization problem in order to find the smoothest frame field interpolating between user prescribed direction constraints. A parameter will allow control over the orthogonality of the field.

Variables. To ease the resolution of the optimization problem, we will not directly use the angle α, β or the vector u, v as variables. Instead, we use the double angle vectors ū = (cos(2α), sin(2α)) and v = (cos(2β), sin(2β)). The vector ū is mapped to either u or -u so the change of variables is valid. With these new variables the coefficients in Eq. ( 4) now read:

c 2,1 (ū, v) = 4ū 1 + 4v 1 c 2,2 (ū, v) = 4ū 2 + 4v 2 c 4,1 (ū, v) = 2ū 2 1 + 2v 2 1 -2 c 4,2 (ū, v) = 2ū 1 ū2 + 2v 1 v2 (5)
Direction constraints. The most common constraint for frame field design is to impose constant frames at certain faces. In our framework, setting these constraints is straightforward : one can compute the angle of each target direction and constraint ū and v to the double angle vectors.

A more interesting constraint, often used for non-orthogonal quad remeshing (see Fig. 5), is to fix only one direction out of two. To do so, one can constrain either ū or v to the desired direction. Note that the choice of ū or v has no impact on the results as the polynomial representation is invariant by reordering of the directions. Constraining either ū or v does not affect the second vector of the frame which can rotate freely.

However, when only one direction is fixed, the second direction can become degenerate if no unit constraint forces it to be non-zero. To avoid such phenomenon, we add, for boundary triangles only, a unit norm energy:

E b (ū, v) = i∈∂F |ū i | 2 -1 2 ,
where ū is the unconstrained direction. This energy is weighted by the parameter b.

For our quad meshing results, we constrain the frames at the boundary to have one direction tangent to the boundary edge. If a triangle has two boundary edges then the frame is fully constrained. In introduction, we emphasis the challenge of remeshing models with sharp corners such as in Figures 4 and11. To catch these geometric features, it is important that the directions are fully constrained where several feature lines meet with acute or obtuse angles.

Smoothness energy. Following Section 1.3, the smoothness energy E s on a 2D frame field depends only on the coefficients of the decomposition Eq. ( 5). We slightly modify the energy to add a parameter λ balancing the frequency of order 2:

E s (ū, v) = λ (i j)∈E m=1,2 c 2,m (ū i , vi ) -c 2,m (ū j , v j ) 2 + (i j)∈E m=1,2 c 4,m (ū i , vi ) -c 4,m (ū j , v j ) 2 .
Orthogonality parameter. As mentioned in Section 2.2, the orthogonality is directly related to second order frequencies. In the smoothness energy, we weight differently the coefficients of order 2, thus putting increasing penalty in these frequencies will lead to frames with constant angle over the entire mesh. If the field have singularities, due the topology of the domain, the frames are forced to be orthogonal. Indeed, a non-orthogonal frame possesses only a π rotation symmetry. So for a singularity to appear, the field must either rotates and change angle at the same time or have index 1/2 singularities. Both these solutions are far more costly in E s than having orthogonal frames. Therefore, introducing a simple weight λ on second order coefficients naturally favors orthogonal frames around singularities and by extension, over the entire shape. The effect of the orthogonality parameter λ on the frame field and the quad-meshes is illustrated in Figure 3. The value λ = 1 correspond to the geometric notion of smoothness explained in Section 1.3. A value greater than one promotes orthogonality and the frame field becomes increasingly similar to the orthogonal frame field. A λ lower the one gives more freedom to the frame field which can be degenerate. As a result, the singularities are pushed on the boundary of the domain.

The choice of this parameter is dictated by the application and is often the result of a trade-off. On one hand, orthogonal frames lead to quad-meshes with elements close to perfect square but can not deal with sharp corners, omnipresent in CAD models. On the other hand, non-orthogonal frame fields compensate sharp features by introducing shearing and on extreme cases removing singularities from the interior of the domain. Optimization. Our optimization scheme is based on the MBO algorithm introduced in [START_REF] Viertel | An approach to quad meshing based on harmonic cross-valued maps and the ginzburg-landau theory[END_REF] to optimize singularity position in cross fields. The idea is to alleviate the unit norm constraint on ū and v by alternating between a diffusion step and renormalization step. Namely, the pair of variables (ū, v) n+1 at iteration n + 1 are obtained by minimizing the smoothness energy while staying close to the iteration n:

(ū, v) n+1 = arg min x,y E s (x, y) + bE b (x, y) + τ∥(x, y) -(ū, v) n ∥ 2 2 , (6)
and then renormalizing (ū, v) n+1 to unit norm. This scheme proves to be the most efficient in term of smoothness of the field and singularity placement. The minimization step is done with LBFGS [START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF].

Results

We assess the usability of our frame field design method in the context of quad-meshing. All the frame fields are constrained to have one direction tangent to the boundary of the domain. To avoid zero norm solutions, the energy promoting unit norm at the boundary is weighted by b = 10. In all our experiments, the diffusion parameter is set to τ = 1/2. The orthogonality parameter λ is tuned according to our need but the value 1.5 proves to be a good compromise in most cases. The optimization problem being non-convex, a good initialization is crucial. We initialize with an orthogonal frame field [START_REF] Ray | N-symmetry direction field design[END_REF] to avoid completely degenerate solutions. The choice of the initial field and its impact on the solution of the optimization are discussed for the three-dimensional case in Section 3.3 but the conclusion also applies in 2D. Our optimization typically needs 10 iterations to converge.

The quad-meshes are extracted from a frame field using the state-of-the-art methods [START_REF] Kälberer | Quadcover-surface parameterization using branched coverings[END_REF][START_REF] Bommes | Integer-grid maps for reliable quad meshing[END_REF].

Quad meshing of CAD models. Generating quad meshing for CAD models is made challenging by the presence of sharp cor- Table 1: Computation time in second for each step of our method: initialization with the orthogonal frame field from [START_REF] Ray | N-symmetry direction field design[END_REF], solving our non-convex optimization problem and the quad mesh extraction from [START_REF] Bommes | Integer-grid maps for reliable quad meshing[END_REF].

ners: at such places a perfect square cannot fit. Figure 4 illustrates this issue. The orthogonal frame field in Fig. 4a is unable to satisfy the constraints on both side of the sharp angle, thus failing to produce a valid quad-mesh. Our method generate anisotropic quad-mesh (Fig. 4b) naturally fitting the sharp corner constraints.

Figure 5 shows a collection of quad-meshes generated from our non-orthogonal frame fields, see how the boundary constraints are nicely satisfied with a small number of singularities while staying away from highly degenerate quads.

Timings. Table 1 shows the timings of each step of our method for each figure of the paper. The most time consuming step is the computation of the non-orthogonal frame field as it requires to solve a non-convex optimization problem. Our algorithm is implemented in C++ but has not been optimized to reach high performances. For instance, we remarked that waiting until convergence of the LBFGS step is often unnecessary and the optimization can be stopped earlier. Also, a Newton method could be used. All the experiments are conducted on a laptop with a four-core, 2.6 GHz Intel Core i5-10400H and 32 GB of RAM.

Comparison with PolyVector field [20]

Diamanti et al. [START_REF] Diamanti | Integrable polyvector fields[END_REF] uses a representation of frame fields as coefficients of a complex polynomial. Like ours, their representation is independent of the direction set ordering. Moreover, they introduce an energy promoting the integrability of the field, i.e. such that each direction is the gradient of a function. Like us, they directly use direction vectors as variables in their numerical scheme. Thus, their integrability energy could be immediately added to our optimization.

Distance between frames. Leaving aside numerical optimization aspects, the main difference resides in the smoothness energy. The representation used in Diamanti et al. [START_REF] Diamanti | Integrable polyvector fields[END_REF] obfuscates the geometrical meaning of the frame representation, leading the authors to seek smoothness of the coefficients of the complex polynomial instead of the smoothness of the field. In our notation, the PolyVector representation of a unit frame boils down to replacing our coefficients c by c:

c2,1 (α, β) = 2 cos(2α) + 2 cos(2β) c2,2 (α, β) = 2 sin(2α) + 2 sin(2β) c4,1 (α, β) = 2 cos(2α + 2β) c4,2 (α, β) = 2 sin(2α + 2β) . (7)
If the second order coefficients matches ours in Eq. ( 4), the fourth order are significantly different. To compare the distance defined by our coefficients c and the coefficients c in Eq. ( 7), we plot the distance between a frame and all its rotation in Figure 6. For orthogonal frames, both distances are identical (Fig. 6a). For non-orthogonal frames the distances differs significantly on two Figure 6: Polar plot of the distance between two orthogonal frames (Fig. 6a) and two non-orthogonal frames (Fig. 6b) with respect to the rotation angle. The distance defined by our representation is plotted in blue and the PolyVector [START_REF] Diamanti | Designing npolyvector fields with complex polynomials[END_REF] distance in red. Insets highlights frames positions at local extrama. The PolyVector distance is not geometrically sound as the maximum distance is reached for overlapping non-orthogonal frames. aspects. First, our distance has globally larger values, thus, with λ = 1 we penalize more the change of angle between frames than PolyVector. Second, the PolyVector distance does not reach its maximum at angle π/2 but at an angle where the frames overlap (Fig. 6b).

Frame field comparison. Using the change of variables (α, β) → (ū, v) in Eq. ( 7):

c2,1 (ū, v) = 2ū 1 + 2v 1 c2,2 (ū, v) = 2ū 2 + 2v 2 c4,1 (ū, v) = 2ū 1 v1 -2ū 2 v2 c4,2 (ū, v) = 2ū 1 v2 + 2ū 2 v1
, we can generate PolyVector fields by replacing c by c in our smoothness energy E s . Figure 7 compares quad-meshes generated from our nonorthogonal direction fields with the one obtained from the coefficients c. We evaluate our method with parameter λ = 0.5 and λ = 1, all other parameters being equal. As expected with λ = 1, our method leads to a quad-mesh with an additional singularity in the upper right fillet that is not present in the PolyVector method [START_REF] Diamanti | Integrable polyvector fields[END_REF]. We found that both methods leads to similar results when we set λ to 0.5. Indeed, although the PolyVector frame distance does not follow a geometric intuition, the difference with our smoothness energy is most prominent when the the frames are faraway from each other. Therefore, it does not greatly impact design of smooth fields as adjacent frames are close to each other.

3D frame fields

In the particular case of 3D frame fields, each frame has three directions noted u, v, w. Instead of the Fourier basis, the space of polynomial P 2k,3 is decomposed with the Spherical Harmonic basis. We will see that when using polynomial of degree 4 our representation of 3D cross field only needs fourth order Spherical Harmonic functions much like the most used representation [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF].

As noted in [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF][START_REF] Ray | Practical 3d frame field generation[END_REF], 3D cross fields can be represented by the spherical polynomial of degree four: p(s) = s 4 1 + s 4 2 + s 4 3 . This is, in fact, the polynomial of lowest degree presenting the same symmetry as a 3D cross. However, higher order polynomials of even degree could be consider for frame field generation. In order to keep the computational complexity as low as possible, we will use degree 4 polynomials.

In this case, the representation polynomial adapted from Eq. ( 1) reads:

p(s) = ⟨u, s⟩ 4 + ⟨v, s⟩ 4 + ⟨w, s⟩ 4 , (8)
and is plotted in two configurations in Figure 8.

Spherical Harmonic decomposition

The Spherical Harmonic functions constitute an orthogonal basis for spherical functions. In computer graphics, they have notably been used for BRDFs approximation for real-time rendering [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF][START_REF] Green | Spherical harmonic lighting: The gritty details[END_REF].

Spherical Harmonics functions Y m ℓ : S 2 → R are identified by two indices : ℓ the band index, corresponding to the frequency of the signal, and m spanning all functions within the frequency ℓ. Any spherical function f : S 2 → R is uniquely decomposed as:

f (s) = ∞ ℓ=0 ℓ m=-ℓ f ℓ,m Y m ℓ (s).
Zonal harmonics. Decomposing the representation polynomial of Eq. ( 8) in the Spherical Harmonics basis is made simple due to their rotational invariance. In particular, the spherical polynomial f : s → s 4 3 , representing the direction z, is invariant by rotation around the z-axis. Thus, its decomposition in SH only uses the functions Y 0 ℓ , ℓ ≥ 0: f (s) = 63Y 0 0 + 36Y 0 2 (s) + 8Y 0 4 (s). Such functions are called zonal harmonic function and have one fundamental property: their rotation to an arbitrary unit axis u is easily computed in the SH basis [START_REF] Sloan | Stupid spherical harmonics (sh) tricks[END_REF]:

⟨u, s⟩ 4 = 63Y 0 0 Y 0 0 + 36 2 m=-2 Y m 2 (u)Y m 2 (s) + 8 4 m=-4
Y m 4 (u)Y m 4 (s). Thus, the decomposition of the frame field polynomial in Eq. ( 8), is simply the sum of the above decomposition:

c 2,m (u, v, w) = Y m 2 (u) + Y m 2 (v) + Y m 2 (w), c 4,m (u, v, w) = Y m 4 (u) + Y m 4 (v) + Y m 4 (w). (9)
Very much like 2D frames, 3D frames only have coefficients in band ℓ = 2 and ℓ = 4. One can check from the expression of the SH that a frame with uniform length (i.e. |u| = |v| = |w|) is orthogonal if and only if all the coefficients c 2,m are zero. So, for orthogonal frames we recover the basic property used when designing 3D cross fields [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF][START_REF] Ray | Practical 3d frame field generation[END_REF]: all rotations of the polynomial i=1,2,3 s 4 i are spanned by the SH of band 4. However, Equation ( 9) is more precise and gives the exact relation between the SH coefficients and the frame deformation.

Designing 3D frame fields with constraints

Unlike the 2D case, there is no obvious change of variables simplifying the expression of the coefficients in Eq. ( 9). Therefore, we keep the three unit vectors u, v, w ∈ S 2 as variables. Like previous works [START_REF] Huang | Boundary aligned smooth 3d crossframe field[END_REF][START_REF] Ray | Practical 3d frame field generation[END_REF][START_REF] Palmer | Algebraic representations for volumetric frame fields[END_REF], performing computation in the space of direction fields requires to solve a non-convex optimization problem.

Direction constraints. Direction constraints are easily set by directly constraining u,v or w. We can fix one, two or three directions at a time. The choice of which of these variables is fixed does not impact the results as our representation is agnostic to reordering of the direction set. This is very convenient for boundary conditions in global parametrization problems.

In all our experiments, a frame associated to a boundary facet has its directions constrained by vectors normal to the boundary. If a tet has one, two or three boundary triangles then the frame must satisfy one, two or three direction constraints.

Smoothness energy. The smoothness energy is a simple distance between coefficients of the decomposition of two adjacent frames as explained in Section 1.3:

E s (u, v, w) = λ (i j)∈E 2 m=-2 c 2,m (u i , v i , w i ) -c 2,m (u j , v j , w j ) 2 + (i j)∈E 4 m=-4 c 4,m (u i , v i , w i ) -c 4,m (u j , v j , w j ) 2 .
Orthogonality parameter. As noted in Section 3.1, the orthogonality of the frame is related to the coefficients of the second band. Thus, we penalize the band coefficients ℓ = 2 in the energy with a weight λ. Increasing λ means that it is harder for the field to change its angle and favor fields with constant angle between its directions. As singularities have smaller energy when the field is orthogonal around singularities, this simple parameter on the order 2 coefficients forces global orthogonality of the field.

Figure 10 illustrates the effect on parameter λ on hexmeshes generation. Values lower than one generates polycube like parametrization with the singularities are at the boundary of the domain while values greater than one approach the cross field solution generated by [START_REF] Ray | Practical 3d frame field generation[END_REF].

Hex-meshing of CAD models

We evaluate the quality of our 3D frame fields by their ability to generate hex-meshes. In all our experiments, hex-meshes are extracted from the frame field using the state-of-the-art method [START_REF] Nieser | Cubecover-parameterization of 3d volumes[END_REF][START_REF] Lyon | Hexex: Robust hexahedral mesh extraction[END_REF]. The orthogonal frames are generated using [START_REF] Ray | Practical 3d frame field generation[END_REF].

The optimization is done the same way as for 2D frame fields. We iterate between minimizing Eq. ( 6) with a diffusion parameters τ = 1, b = 100 and renormalization of the vectors u, v and w.

The presence of sharp corners is a classic challenge for hexremeshing methods. This situation occurs frequently for CAD model and is often badly handled by cross fields which collapses the hexes at the corner as in Figure 11a. Whereas the direction field designed by our method introduces a singularity and a valid hexmesh can be extracted (Figure 11b). To produce such results, it is primordial that frames at sharp corners are fully constrained.

Figure 9 also highlights the necessity of non-orthogonal frame fields to capture sharp corners. Cross fields generates hexmeshes with heavily distorted or missing elements while our method is able to fit a valid hex element in the corner.

In Figure 12, the presence of shear forces the orthogonal frames to oscillate to satisfy the boundary constraints while our frame field can remain constant.

Initialization. The energy E s is non-convex, thus its outcome is highly dependent of the initialization. Figure 13 compares solutions of our optimization problem for three different starting points. The initialization with a random orthonormal frame field always converges to the same solution (Fig. 13b). An initialization with a smooth cross field [START_REF] Ray | Practical 3d frame field generation[END_REF] needs less iterations to reach convergence but the solution is still the same (Fig. 13a). As reported in Table 2, this reduces (sometimes considerably) the computation time. However, our optimization scheme gets trapped in a local minimum when initiated with a random field with nearly flat frames (Fig. 13c). At convergence, the total energy is higher than the optimal solution. Moreover, the field is corrupted by degenerate frames (circled in red) making it useless for meshing applications.

In all our experiments, we initialize the optimization with the smooth orthonormal frame field obtained by [START_REF] Ray | Practical 3d frame field generation[END_REF]. Table 2: Computation time in second for each step of our 3D frame field design method: initialization with an orthogonal field [START_REF] Ray | Practical 3d frame field generation[END_REF], solving our non-convex optimization problem and the hex-mesh extraction using [START_REF] Lyon | Hexex: Robust hexahedral mesh extraction[END_REF].

Timings. Table 2 lists the computation time for generating a nonorthogonal frame field with our method. Again, the most time consuming step is the computation of the non-orthogonal frame field as it requires to solve a non-convex optimization problem.

All the experiments are conducted on a laptop with a four-core, 2.6 GHz Intel Core i5-10400H and 32 GB of RAM. All algorithms are implemented in C++.

Conclusion

In this paper, we introduced a representation of direction fields common to surfaces and volumes, allowing constraints on one, two or three directions independently with an intuitive control on the frame orthogonality. We demonstrated that it can be useful for many geometry processing applications such as anisotropic planar and volumetric remeshing.

While our optimization is non-convex and more challenging than standard cross field representation, it allows more flexible constraints which are better suited for global parametrization problems.

In the future, we plan to study integrable direction fields, i.e. fields whose directions are gradients of scalar functions. This is fundamental to guarantee the existence of a bijective global parametrization. This work has been conducted for fields on surfaces but is still missing for volumetric meshes. With our representation we can hope to extend two dimensional results to 3D. The initialization with a completely random field, including degenerate frames, gets trapped in a degenerate local minimum, resulting in a higher total energy E (Fig. 13c).

Figure 1 :

 1 Figure 1: A unique direction (black arrows) represented by a polar (left) and spherical (right) polynomial of degree 2. The value of the polynomial at each point of the circle is represented as a variation of the radius in polar coordinates.

  (a) Orthogonal frame (b) Non-orthogonal deformation

Figure 2 :

 2 Figure 2: A 2-Direction set {u, -u, v, -v} (black arrows) represented by a polar polynomial of degree four.

  (a) λ = 0.1 (b) λ = 0.5 (c) λ = 1 (d) λ = 1.5 (e) Orthogonal [25]

Figure 3 :

 3 Figure 3: Results of our frame field design method for increasing values of the orthogonality parameter λ. Low values (left) push the singularities out of the domain while high values approaches the cross field obtained by Viertel et al. [25] (right).

Figure 4 :

 4 Figure 4: An orthogonal frame field (4a) cannot be used to generate a quad-mesh for a model with sharp corners. Our non-orthogonal frame field (4b) provides a clean solution.

Figure 5 :

 5 Figure 5: Our method is able generate nicely shaped quad-meshes by adjusting the parameter λ = 1.5.

  Fig. FacesInit Opti Meshing3 658 0.017 <1 2 4b 1282 0.026 <1 3 5c 9968 0.177 26 13 5a 18430 0.355 66 30 5d 29000 0.611 120 45 5b 36047 0.775 161 46

1 Figure 7 :

 17 Figure 7: The smoothness energy from PolyVector field [20] achieves similar results than ours for λ = 0.5.

Figure 8 : 3 -

 83 Figure 8: 3-Direction sets {u, -u, v, -v, w, -w} represented by spherical polynomials of degree 4.

Figure 9 :

 9 Figure 9: Extruded square with curl. The first result uses an orthogonal field, the second uses our fields. The following two images simply add more curl to the extrusion. Sharp corners (disk highlights) are better captured with our fields. The layer of hexahedra inside the volume also remains parallel to the extremities.

Fig

  Fig. Tets Init Opti Meshing 13a 862 0.016 22 0.7 13b 862 -25 0.7 13c 862 -26 -9 4420 0.096 31 3.0 10 13225 0.349 454 8.4 12b 10788 0.332 305 7.7

  (a) λ = 0.7 (b) λ = 1 (c) λ = 5 (d) Orthogonal [21]

Figure 10 :

 10 Figure 10: Results of our 3D frame field design method for increasing values of the orthogonality parameter λ. Low values (left) push the singularities out of the domain while high values approaches the 3D cross field parametrization obtained by Ray et al. [21] (right).

1 Figure 11 :

 111 Figure 11: On a CAD model with a sharp corner, the 3D cross field collapses part of the model (Fig. (11a)). Our method successfully adds a singularity to compensate the slope (Fig. (11b)).

  (a) Cross field[START_REF] Ray | Practical 3d frame field generation[END_REF] (b) Ours λ = 1

Figure 12 :

 12 Figure 12: On a CAD model with heavy shear, the 3D cross field from creates unnecessary oscillation to compensate (Fig. (12a)). Our method outputs nicely shaped and regular hexes (Fig. (12b)).

Figure 13 :

 13 Figure13: Results of our 3D frame field design method starting from three different initializations. Starting from a smooth orthogonal frame field[START_REF] Ray | Practical 3d frame field generation[END_REF] (Fig.13atop) or from a random orthogonal field (Fig.13btop) leads to the same solution (Fig.13a-Fig.13b bottom). The initialization with a completely random field, including degenerate frames, gets trapped in a degenerate local minimum, resulting in a higher total energy E (Fig.13c).

Preprint submitted to ElsevierJuly 15, 2021