Maria Maria
email: hmaria@physik.uzh.ch

Han Veiga

Remi Abgrall
email: remi.abgrall@uzh.ch

Maria Han Veiga

Rémi Abgrall

Towards a general stabilisation method for conservation laws using a multilayer Perceptron neural network: 1D scalar and system of equations

Keywords: limiters, computing methods, neural networks

published or not. The documents may come

Towards a general stabilisation method for conservation laws using a multilayer Perceptron neural network: 1D scalar and system of equations

Introduction

When dealing with nonlinear conservation laws [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], it is well known that discontinuous solutions can emerge, even for smooth initial data. The numerical approximation of the discontinuous solution will develop nonphysical oscillations around the discontinuity, which in turn will negatively impact the accuracy of the numerical scheme. There are different methods to control these oscillations, namely, adding a viscous term as in [START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] or limiting the solution.

∂ ∂t u + ∇ • f (u) = 0 (1)
∂ ∂t u + ∇ • f (u) = ∇ • (ν(u)∇u) (2)
Neural networks gained new popularity due to the computational tractability of the back-propagation algorithm, used for the learning of weights and biases in a deep neural network. Furthermore, it has been empirically shown to generate robust models for classification in many areas of applications [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Sutskever | Sequence to sequence learning with neural networks[END_REF] and theoretically, to generate universal classifiers and universal function approximators [START_REF] Rojas | Deepest Neural Networks[END_REF][START_REF] Rojas | Networks of width one are universal classifiers[END_REF][START_REF] Sutskever | Deep, narrow sigmoid belief networks are universal approximators[END_REF].

We can write a stabilisation method as a function f that takes some local properties of the solution (denoted by the map X(u(x))), and returns a modified solution value ũ(x), which has some desired properties, e.g. non-oscillatory, maximum principle preserving. f u(x), X(u(x)) = ũ(x)

For example, in the context of finite volume schemes, we can write the Minmod limiter as: f ((u j ), X(u j )) = u j -minmod(X(u j )) = ũj where the map X(u j ) encodes local properties of the solution at u(x j ). In this case, X(u j ) = (ū ju + j-1/2 , ∆ + ūj , ∆ -ūj ), where the usual notation applies: ū the cell average, ∆ + ūj := u j+1 -u j , ∆ -ūj := u j -u j-1 and u * = u(x * ).

It is well known this clips the extrema for smooth solutions. To minimize this effect, one can consider a TVD limiter [START_REF] Arora | A well-behaved tvd limiter for high-resolution calculations of unsteady flow[END_REF]:

f (u(x j ), X(u(x j )), M ) = u(x j ) -TVD(X(u(x j )), M ) = ũj where the map X(u(x j )) = (ū j -u + j-1/2 , ∆ + ūj , ∆ -ūj ) again encodes local properties of the solution at u(x j ), and 2 , where M is a user defined parameter that gives an estimate of the smoothness of the solution u(x). It is to note that M can take a value in a large range of positive numbers, and that it is usually a global quantity fixed in the beginning of the numerical experiment. Thus, this can be a drawback if the solution has different smoothness properties across the domain.

T V D(X(u j ), M ) = minmod(X(u j )) else u(x j ) if |ū j -u + j-1/2 | ≤ M (∆x)
In the context of discontinuous Galerkin method, we can mention, among many, the high-order limiter by [START_REF] Krivodonova | Limiters for high-order discontinuous Galerkin methods[END_REF], which does the limiting in a hierarchical manner, relying on the modal representation of the solution u(x). This does not clip the solution maxima, thus having desirable properties, but it is formulated specifically for the modal discontinuous Galerkin method.

In this work we learn a data-driven stabilisation function f (X(u(x))) which doesn't require user input, and we explicitly show how to integrate such stabilisation method with a working 1-dimensional computational fluid dynamics (CFD) code. Furthermore, we design this stabilisation function with the intent to be independent of the underlying numerical scheme, so the input required is formulated by nodal values of the solution and finite differences.

The paper is structured as follows: in section 2 the construction of the dataset is described, the methodology in presented in section 3, numerical results are shown in section 4 and conclusions in chapter 5.

In the spirit of open science, all the codes and trained models are made available in a public repository [8].

Dataset

The dataset is an integral part of data-driven analysis. It contains the data for which we want to learn a mapping for. For our task, to learn a function which indicates whether a discontinuity is present in the solution or not, our data is the set containing N samples {X i , y i } N i=1 , where X i denotes some local properties of the solution u i (features) and y i (labels) indicates the existence (or not) of a discontinuity.

The dataset is generated by performing many runs on a 1-dimensional discontinuous Galerkin code solving the advection equation for different initial conditions, orders and mesh sizes (see table 1), and the labels are obtained by running a high-order limiter for discontinuous Galerkin [START_REF] Krivodonova | Limiters for high-order discontinuous Galerkin methods[END_REF] which does hierarchical limiting of the modes. If the limited solution deviates from the unlimited solution more than ζ, then the cell needs stabilisation, and we get a positive label. We set ζ = 0.01.

The interesting point of generating the labels in this way is that we can learn a black-box function of a limiter which is specific to a particular framework (in this case, discontinuous Galerkin) and integrate this with other numerical methods. Furthermore, the cell which are flagged by the high-order limiter are similar to a well tuned TVD limiter for approximation orders higher than 2.

Features

The features X i are the different quantities used to describe the local solution u i . For the sake of generalisation, we choose features which are readily available in different numerical methods, such as: cell mean value, values at interface, divided differences between neighbours, etc. The complete description of the features can be found in table 2.

Of course, one could consider features which are specific for a particular method, e.g. the modes of the solution when using a discontinuous Galerkin method, or the smoothness parameters for a WENO method.

ID Feature Name

Method

In the following section, we focus on the three main aspects of this method:

1. the set-up of the learning algorithm;

2. the integration of a trained model with an existing CFD code;

3. the performance measures used to validate this method.

For the specific implementation details, please refer to the github repository [8].

Training a neural network offline

In this section we describe the details of the learning algorithm. We use a multilayer Perceptron neural network. For a general comprehensive introduction of neural networks one can refer to [START_REF] Goodfellow | Deep Learning[END_REF]. We wish to learn a map f : X → Y, where X denotes an arbitrary set containing examples that we wish to label with possible outcomes Y. The task at hand is a binary classification (i.e. does this cell require stabilisation?), thus f will be our binary classifier and Y = {0, 1}.

We choose f to be defined by the composition of a sequence of functions g 1 , g 2 , ...g n , yielding

f (x) = g n (...g 2 (g 1 (x))).
This is known as a deep neural network. There are many different classifiers which can be used, but it has been shown that deep neural networks perform well on a variety of classification tasks, in particular when the classification plane is nonlinear [START_REF] Goodfellow | Deep Learning[END_REF].

Each function

g i (w i , b i , h i (•)
) is parametrized by a matrix w i , called weights, a vector b i called bias and an activation function h i (•) which introduces the non-linearity on the neural network.

These parameters are tuned through a loss function L(x), which measures how well the mapping f performs on a given dataset D, using the gradient descent algorithm and back-propagation.

The gradient descent [START_REF] Snyman | Practical Mathematical Optimization: An Introduction to Basic Optimization Theory and Classical and New Gradient-Based Algorithms[END_REF] is a first-order iterative optimization algorithm for finding the minimum of a function (in this case, the loss function L(x)), relying on the fact that for small enough η:

a n+1 = a n -η∇L(a n ) then L(a n+1 ) ≤ L(a n ),
for a defined and differentiable loss function L. The size of the update on the opposite direction of the steepest gradient is tuned through η, the learning rate. Instead of using a global learning rate, we use the Adam algorithm [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] to adaptively estimate the update for each parameter (in this case, for the weights and biases), which uses information from the first and second moments of the gradient ∇L(a n ) (with respect to parameter a). The update follows:

a n+1 = a n -α mn √ vn + ε ,
where α denotes the stepsize, typically set to α = 0.001, mt denotes the corrected moving first moment of the gradient, vt the corrected moving second moment of the gradient (both with respect to parameter a) and ε a small value to avoid division by zero. For the full algorithm, please refer to appendix B. Finally, the back-propagation algorithm is used, which gives a rule to update the different w i and b i in order to minimize the loss function [START_REF] Feldman | Neural Networks: A Systematic Introduction[END_REF].

Finally, two different loss functions are used:

• cross-entropy:

L(D) = - 1 N N i y i log(p i ) + (1 -y i ) log(1 -pi ). (3)
• weighted cross-entropy:

L(D) = - 1 N N i y i log(p i )ω + (1 -y i ) log(1 -pi ). (4)
One particular property of our problem is that it is expected that there will be a class imbalance on the dataset (both during the training phase and prediction phase). In particular, it is more likely to find cells which are in no need for stabilisation than ones which are in need for stabilisation. Furthermore, it is more desirable to overlimit than to miss a cell that needs limiting, as it might lead to unphysical results and potentially crash the code. To bias the learning towards predicting the less represented class (needing stabilisation), it is common practice to use a weighted cost function, which increases the penalty of mislabeling a positive label [START_REF] Hoens | Imbalanced Datasets: From Sampling to Classifiers[END_REF]. The asymmetry in the loss function is added through the coefficient ω.

Finally, we specify the activation functions used: for the hidden layers 1, ..., n -1 a Rectified Linear Unit (ReLU) is used:

h(x) = max(0, x).
Although there exist more sophisticated activation functions, typically modifications to ReLU, e.g. leaky ReLU, Parametric ReLU or Randomised leaky ReLU, these require further parameter estimation, add at least one more dependence to the saved model, and the empirical improvement on the performance is not extremely significant [START_REF] Xu | Empirical evaluation of rectified activations in convolutional network[END_REF]. Furthermore, there are theoretical predictions when using ReLU as an activation function, which we make use in the following section to roughly estimate the size of the neural network.

For the last layer (output layer), a sigmoid function is used:

h(x) = 1 1 + exp -x ,
in order to attain a value that can be interpreted as a probability.

Architecture

In this work we trained a fully connected, multi-layer Perceptron (MLP). We follow the results in [START_REF] Petersen | Optimal approximation of piecewise smooth functions using deep ReLU neural networks[END_REF], which establish lower bounds for the number of required non-zero weights and layers necessary for a MLP using ReLU activation functions used, to approximate the classification maps

f : [-1/2, 1/2] d → Y: f = k≤N c k X K k ,
where f is a piecewise smooth function with N "pieces". Here each X K k denotes an indicator function:

X K k (x) = 1, if x ∈ K k 0, if x ∈ K k (x, t) ∈ [0, 1] × R + , (5)
for a compact set

K k ⊂ [-1/2, 1/2] d with ∂K k ∈ C β , the set of β times differentiable curves, for all k ≤ N .
Namely, β controls the regularity of the boundary of the set K k , d denotes the dimension of the feature map. Theorem 4.2 [START_REF] Petersen | Optimal approximation of piecewise smooth functions using deep ReLU neural networks[END_REF] provides the following lower bound for the minimal number of nonzero weights that a neural network R ρ (Φ f ε ) needs to have in order to approximate any function belonging to a function class C ∈ L 2 ([-1/2, 1/2] d ), up to an L 2 error of at most ε:

#W eights ≥ C • ε -2(d-1)/β / log 2 1 ε , where C = C(d, β, B, C 0 ) is a constant, B, C 0 > 0,
B is an upper bound on the supremum norm of R ρ (Φ f ε ) and C 0 is related to with the maximum bits needed to encode a weight of the network.

And Theorem 3.5 [START_REF] Petersen | Optimal approximation of piecewise smooth functions using deep ReLU neural networks[END_REF] provides approximation rates for piecewise constant functions. In particular, the number of required layers depends on the dimension d and for the regularity parameter β:

#Layers ≥ c log 2 (2 + β) • (1 + β/d).
Making optimistic assumptions about the classification function and the desired approximation accuracy (in L 2 norm), we can bound, from below, the number of necessary neurons and layers as (non-zero) #W eights ≈ 2.7 × 10 8 and #Layers ≈ 2.4, respectively, for d = 11, β = 2 and ε = 0.10. As these are lower bounds and C, c are unspecified, we fix the number of neurons to be at least 2 28 which is distributed across the varying number of layers. In appendix A, the explicit architectures are detailed.

Integration of the method on a CFD code

We assume that exists a neural network which has been trained offline. Each layer can be fully characterized by the following information:

g i (w i , b i , h i (•))
where h i denotes the activation function, w i the weights matrix and b i a bias vector. Thus, there are two routines necessary to integrate a trained neural network with an existing code:

1. Generation of features: given the local solution u, generate the feature quantities X(u) as specified on section 2.

2. Prediction routine: Given the features X(u), it is necessary to evaluate the function f . As denoted previously, a layer g i is fully specified by parameters w i , b i , h i (•). Once the neural network has been trained offline, the weights and biases can be loaded onto the CFD code. What remains to be implemented are the activation functions for the hidden layers and the activation function for the output layer in order to evaluate f at some given input.

Then, the stabilisation algorithm reads: Data: solution at cell i, u i Result: stabilised solution ũi X = generateFeatures(u i ); label = predict(X); if label = 1 then ũi = stabilisation(u i ); else ũi = u i ; end Algorithm 1: Limiting For systems, we tried both limiting on conservative variables and primitive variables. Each variable is limited independently.

Measuring performance of model

In this work we use two sets of performance measures, namely:

1. Label prediction measures 2. L 1 -norm of the solution For the first set of measures, we can use typical metrics used in computational statistics and machine learning communities:

tp + tn tp + tn + f p + f n (accuracy) tp tp + f n (recall) tp tp + f p (precision)
where tp is the number of correctly predicted positive labels, tn the number of correctly predicted negative labels, f p the number of incorrectly predicted positive labels and f n the number of incorrectly predicted negative labels.

It is important to consider recall and precision because the distribution of labels can be imbalanced. If only accuracy is measured, it is possible to attain a high accuracy by solely predicting the majority label.

For the second set of measures, we consider the L 1 -norm because it is the relevant measure for a CFD code and we can study the effect of this method on the empirical error of the numerical solution.

Implementation details

The neural network was set up and trained using Tensorflow with a python API. The numerical solver is written in Fortran, based on [ a θ denotes the estimate of the ratio between positive and negative labels, which could be measured empirically through the training data. Empirically, we have θ = 0.1

Numerical experiments

This section is split in two parts: in subsection 4.1, we show the performance of several trained neural networks on an unseen validation set by measuring the accuracy, recall and precision. The second part (subsections 4.2 and 4.3), we choose a model that performed well and we integrate it with a CFD code. The model runs as a blackbox limiter (denoted as NN) and we compare its performance to the Minmod limiter and hierarchical high order limiter (denoted as HIO) through the L 1 error norm. We also compare the performance of this limiter with the optimally tuned TVD limiter, although we do not present convergence rates of the TVD limiter due to the dependency on the M parameter. We perform some tests for the linear advection equation and Euler system of equations. The initial conditions are chosen as different from the ones used for the training.

Detection rate

We measure the performance of several models on a unseen validation set. We measure whether the models are able to predict the right labels which are assigned from hand tuned limiters.

On an unseen validation set, a subset of the dataset generated as describes in section 2.

Going forward, we select model 5 as it performs well and the resulting size of the weights matrices per layer is significantly smaller than model 3. It is debatable whether the differences between models 3, 4 and 5 are statistically significant.

Linear Advection

Consider a linear advection equation with a ∈ R:

∂ ∂t u + a ∂ ∂x u = 0 (6)
and periodic boundary conditions.

Gaussian pulse

We consider the following initial condition:

u 0 (x) = 1 + 3 exp(-100(x -0.5) 2 ) (x, t) ∈ [0, 1] × R + (7)
with advection velocity a = 1.

The convergence is shown in tables 3 and 4 after one full crossing for orders 2 and 3. Furthermore, in figure 1 we show how the maxima is clipped using different methods for a grid size of N = 40. We note that the Minmod limiter clips the maximum value of the solution, whereas a well tuned TVD limiter does not. The hierarchical high order limiter (denoted as HIO) throughout this section behaves as Minmod for the second order case, but for third order it does not limit the solution. The learned limiter (denoted as NN) behaves similarly for both orders, but it is to note that it seems to outperform the HIO limiter for the 4: L 1 error for one crossing of the Gaussian pulse 7 using different limiters for order 3.

Smooth and square hat

Next, we consider the following initial conditions, which contain a smooth Gaussian pulse and a hat function:

u 0 (x) = 2, |x -0.7| ≤ 0.1 1 + exp -(x-0.25) 2 2×0.05 2 , otherwise (x, t) ∈ [0, 1] × R + , (8)
again with advection velocity a = 1.

Euler equation

Now we consider the 1-dimensional Euler equations, which describe the behavior of an inviscid flow. This system of equations describe the evolution of a density ρ, a velocity u, a pressure p and total energy E. 6: L 1 error for one crossing of the Gaussian and hat pulses 8 using different limiters for order 3.

∂ ∂t ρ + ∂ ∂x (ρu) = 0 (9)
∂ ∂t (ρu) + ∂ ∂x (ρu 2 + p) = 0 (10)
∂ ∂t E + ∂ ∂x ((E + p)u) = 0. ( 11
)
The system is closed equation of state for an ideal gas:

ρe = p/(γ -1)
,

where e = E -1 2 ρu 2 is the internal energy.

Sod shock tube

We consider the standard Sod shock tube test, given by the initial conditions: (ρ, v, p)(x, 0) = (1.0, 0.0, 1.0) 0.0 < x ≤ 0.5 (0.125, 0.0, 0.1) 0.5 < x < 1.0

(x, t) ∈ [0, 1] × [0.0, 0.1], (12)
and γ = 1.4 and gradient free boundary conditions.

In figure 3 we show the comparison between different limiters at T = 0.1. We show the performance of a tuned TVD limiter and the NN limiter. We note that for the density field (left panel in figure 4.3.1), the solution produced by the NN limiter seems oscillation free and sharper than the solution by the TVD limiter, but for the velocity field (right panel in figure 4.3.1) some oscillations were not corrected enough, and clearly the behaviour of the TVD limiter is more desirable. This could be due to the fact that the shock detection is done over the conserved variables. This could be improved by performing the shock detection over the characteristic variables.

Blast wave

Next we consider the interacting blast waves test, given by the initial conditions: (ρ, v, p)(x, 0) =    (1.0, 0.0, 1000.0) 0.0 < x ≤ 0.1 (1.0, 0.0, 0.01) 0.1 < x ≤ 0.9 (1.0, 0.0, 100.0) 0.9 < x < 1.0

(x, t) ∈ [0, 1] × [0.0, 0.038] (13)
with γ = 1.4 and reflexive boundary conditions. In figures 4 and 5 we show the comparison between different limiters at T = 0.038 for different orders. The unlimited solution is not shown because the code crashes due to the pressure becoming negative shortly after the start of the simulation for orders higher than 1. The dashed line denotes a high resolution solution, run with N = 1000, third order with the HIO limiter. We can note that the NN limiter is not as good at suppressing oscillations as Minmod and the HIO limiter, but stabilises the solution enough to finish the run. Furthermore, we note that the peak is better preserved, which means that it looks like the limiting is less strong than Minmod and the HIO limiter.

Conclusion

The purpose of this work was primarily to demonstrate the potential of using learning algorithms in CFD codes. In particular, we show the different stages necessary to train a black box limiter that can be integrated in different codes. To this end, we described how to construct a dataset for the training phase, how to train a multi-layer Perceptron that detects shocks and how to integrate the trained model with a 1-dimensional discontinuous Galerkin code to serve as a limiter. Then, we performed the usual validation of the limiter in the context of scalar and systems of equations.

An idea that we tested in this paper was training an algorithm that mimics a limiter designed for a particular numerical scheme. In particular, we used the cells flagged by the hierarchical high order limiter [START_REF] Krivodonova | Limiters for high-order discontinuous Galerkin methods[END_REF] for discontinuous Galerkin method, which relies on the fact that the high order coefficients of the Legendre basis functions serve as proxies to the high order derivatives of the solution to perform the limiting, and we learn the mapping which using only nodal values of the solution as described in table 2. This is hopefully a step towards the development of algorithms which could be numerical scheme independent.

We showed empirically that in principle it's possible to learn a parameter free limiter (after the training phase) from data. While the performance for the seen functions is relatively good (as observed in 4.1), the on-the-fly performance could be improved with the careful design of the loss function during the training phase (for example, to include information on the maximum preserving principle). For the smooth case, we observed that this limiter was by far less diffuse than the minmod limiter. Certain oscillations were corrected but there were other oscillations which were not stabilized enough. When applied to a system of equations, we observed that this method still produced sensible results, slightly better than the unlimited version.

It is to note that the training set used was extremely reduced -we only used data from solving the linear advection equation with different initial conditions and varying mesh sizes. Even so, the direct application of this model to unseen initial conditions, to the Euler system and to differing grid sizes was somewhat successful. In a realistic setting, one would use a much more complete dataset.

Of course one must mention caveats in this type of work. For once, there is overhead when using the neural network limiter, as one can observe in table 7, where the average wall-clock time per timestep using with different limiters is shown. This can be improved significantly by using an efficient matrix multiplication routine to evaluate the prediction for each cell. A more concerning problem is the lack of theoretical guarantees for this type of black-box limiter. For example, additional work must be done to build limiters which are maximum principle preserving. These properties can be added by designing an appropriate loss function. In addition, work must be done for the extension to multi-dimensional problems. Table 7: Average time per iteration (s) for sine advection case While this work is a proof-of-concept, it is our belief that these ideas can be applied to other problems which depend on certain local properties of the numerical solution, ultimately contributing towards Computational Fluid Dynamics (CFD) codes which are robust to different initial conditions and that require less parameter tuning to produce readily usable results.

A Appendix: MLP architectures

Here we detail the architectures tested in this work. For this work we found the lower bound (2 28 ) for the number of weights and fixed the activation functions to be ReLU (as detailed in section 3.1.1). What we observed was that this quantity was not producing very good models. Thus, we used (2 32 ) non-zero weights, which presupposes a constant C of order 10. It is left to specify the distribution of the weights across the different layers. The number of weights and neurons is related by multiplying d through the number of neurons per layers to get the total number of weights.

  010765 | 1.7 0.392397 | 0.2 0.224807 | 0.3 0.354754 | 0.2 20 0.003666 | 2.8 0.143413 | 1.5 0.007157 | 5.0 0.107731 | 1.7 40 0.000192 | 2.2 0.038610 | 1.7 0.002575 | 3.2 0.017128 | 2.2 100 0.000078 | 2.6 0.004917 | 1.9 0.000781 | 2.3 0.001750 | 2.3 Table

  164609 | 0.4 0.333485 | 0.2 0.333485 | 0.2 0.297200 | 0.3 20 0.100032 | 0.7 0.170726 | 1.0 0.170726 | 1.0 0.129082 | 1.2 40 0.053035 | 0.8 0.062677 | 1.2 0.062677 | 1.2 0.065400 | 1.1 80 0.029132 | 0.8 0.028890 | 1.2 0.028890 | 1.2 0.032303 | 1.1 100 0.024246 | 0.8 0.023000 | 1.2 0.023000 | 1.2 0.026439 | 1.0 Table 5: L 1 error for one crossing of the Gaussian and hat pulses 8 using different limiters for order 2. 054616 | 0.6 0.339167 | 0.2 0.307901 | 0.3 0.302808 | 0.3 20 0.032699 | 0.7 0.172127 | 1.0 0.075238 | 2.0 0.130308 | 1.2 40 0.022798 | 0.6 0.063178 | 1.2 0.033082 | 1.6 0.065783 | 1.1 80 0.013853 | 0.6 0.029148 | 1.2 0.017352 | 1.4 0.032443 | 1.1 100 0.011725 | 0.7 0.023202 | 1.2 0.014379 | 1.3 0.026401 | 1.0 Table

Figure 2 :

 2 Figure 2: Maxima clipping of the Gaussian and hat pulses 8 after one full crossing for approximation order of 2 (left) and 3 (right) for grid size N = 40.

Figure 3 :

 3 Figure 3: Density and velocity fields for the Sod shock tube 12 at T = 0.1 for grid size N = 64 and approximation order 2.

Figure 4 :Figure 5 :

 45 Figure 4: Density, velocity and pressure fields of the blast wave interaction 13 at T = 0.038 for grid size N = 100 and approximation order 2.

15 ×

 15 10 -5 8.61 × 10 -5 1.29 × 10 -4 TVD limiter 8.70 × 10 -5 1.17 × 10 -4 1.74 × 10 -4 NN limiter 8.22 × 10 -3 1.68 × 10 -2 3.20 × 10 -2

Table 1 :

 1 Runs used to generate the dataset.

	Initial condition sin(2πx),	4 0.25 ≤ x ≤ 0.75 1 otherwise	,	1 2 sin(2πx) x ≤ 0.3 0 otherwise	x ∈ [0, 1]
	Mesh size	8, 16, 32, 64, 128
	Order	2, 3

Table 2 :

 2 Features table

			Description
	1	h	cell width
	2	ūi	average value of solution at cell i
	3	ūi+1	average value of solution at cell i + 1
	4	ūi-1	average value of solution at cell i -1
	5 6 7 8	u + i-1 2 u -i+ 1 2 u -i-1 2 u + 2 i+ 1	value of solution at interface i -1/2 reconstructed in cell i value of solution at interface i + 1/2 reconstructed in cell i value of solution at interface i -1/2 reconstructed at cell i -1 value of solution at interface i + 1 reconstructed at cell i + 1
	9	du i+1	divided difference between ūi and ūi+1 , divided by h
	10	du i-1	divided difference between ūi and ūi-1 , divided by h
	11	du i	divided difference between ūi+1 and ūi-1 , divided by 2h

  16].

	Model	Description	Accuracy (%) Recall (%) Precision (%)
	Random	randomly guessing	50.00	θ a	50.00
	Model 1	2 hidden layers (HL)	77.96	69.17	24.13
	Model 2	3 HL	95.52	80.66	84.19
	Model 3	4 HL	95.15	79.26	82.67
	Model 4	5 HL	95.14	81.71	81.01
	Model 5 5 HL + weighted loss (ω = 5)	95.67	79.05	86.39

Table 3 :

 3 It makes sense that the performance does not depend (as much as the HIO limiter) on the order of the method, as we train the algorithm with only nodal information. | 0.4 0.389622 | 0.2 0.389622 | 0.2 0.353091 | 0.2 20 0.055585 | 1.8 0.142990 | 1.4 0.142990 | 1.4 0.107405 | 1.7 40 0.008065 | 2.3 0.038609 | 1.7 0.038609 | 1.7 0.017078 | 2.2 100 0.001181 | 2.3 0.005105 | 1.9 0.005105 | 1.9 0.002031 | 2.3 L 1 error for one crossing of the Gaussian pulse 7 using different limiters for order 2.

	3.4 3.6 3.8 4.0					T=0 TVD M=100 NN Unlimited Minmod HIO	3.4 3.6 3.8 4.0					TVD M=100 NN Unlimited Minmod HIO T=0
	3.2							3.2
	3.0							3.0
	2.8							2.8
	2.6							2.6
	0.40	0.45	0.50	0.55	0.60	0.65	0.70	0.40	0.45	0.50	0.55	0.60	0.65	0.70
	Figure 1: Maxima clipping of the Gaussian pulse 7 after one full crossing for approximation order of 2 (left)
	3 (right) for grid size N = 40.
	second order case. N	No Limiter	MinMod		HIO			NN
		10 0.197729

B Appendix: Adam algorithm

The Adam algorithm [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] is used for stochastic optimisation. The algorithm updates exponential moving averages of the gradient (m t ) and the squared gradients (v t ) where the hyper-parameters β 1 , β 2 ∈ [0, 1) control the exponential decay of these moving averages. These moving averages are estimates of the first and second moment (uncentered variance) of the gradient. Because the quantities m 0 , v 0 are initialised at 0, they are biased towards zero (in particular if the decaying rates β 1 and β 2 are small). Hence, there is an additional correction leading to the unbiased quantities mt and vt . The algorithm is detailed below.

Data: α: stepsize, β 1 , β 2 ∈ [0, 1]: exponential decay rates for the moment estimates, f (θ): stochastic objective function with parameters θ, θ 0 : initial parameter vector Result: θ t (resulting parameters) m 0 = 0 (initialise 1st moment vector); v 0 = 0 (initialise 2nd moment vector); t = 0 (initialise timestep); while θ t not converged do t = t + 1;


David González

Icíar Alfaro

Carlos Quesada
email: cquesada@unizar.es

Elías Cueto
email: ecueto@unizar.es

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Computational vademecums for the real-time simulation of haptic collision between nonlinear solids

Keywords: Real time, contact, Model reduction, Proper Generalized Decomposition, Parametric models, Computational vademecums

In this paper a novel strategy is presented for the real-time simulation of contact between non-linear deformable solids at haptic feedback rates. The proposed method is somehow related with the Voxmap Pointshell method for two deformable solids. Its novelty and crucial advantages over existing implementations of this algorithm come from the intensive use of computational vademecums. These are in essence a pre-computed solution of a parametric model in which every possible situation during the on-line phase of the method has been considered through the introduction of the appropriate parameters. Such a high-dimensional parametric model is efficiently solved by using Proper Generalized Decompositions (PGD) and stored in memory as

Introduction

Computational contact mechanics [START_REF] Wriggers | Computational Contact Mechanics[END_REF] constitutes nowadays a very active field of research due to its inherent difficulty, associated to the highly nonlinear character of its models. But when we deal with non-linear solids and, in addition, real-time response is required, the problem becomes extremely burdensome and has generated a plethora of publications searching for a good compromise between accuracy and time to response, see for instance [START_REF] Barbič | Time-critical distributed contact for 6-DoF haptic rendering of adaptively sampled reduced deformable models[END_REF] [START_REF] James | Bd-tree: Output-sensitive collision detection for reduced deformable models[END_REF] [10] [START_REF] Teschner | Collision detection for deformable objects[END_REF], to name but a few of the available references.

Haptic peripherals have become very popular for augmented simulation in immersive environments, particularly for education purposes and games.

They have been used notably in medical environments as an essential tool for the education in minimally invasive procedures [START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF] [60] [START_REF] Delingette | Hepatic surgery simulation[END_REF] [40] [START_REF] Meier | Realtime deformable models for surgery simulation: a survey[END_REF] [18] [START_REF] Alfaro | Real-time in silico exper-iments on gene regulatory networks and surgery simulation on handheld devices[END_REF] [49] and for training of complex industrial processes. In particular, aircraft and automobile industries have recently incorporated virtual reality to the design process to see whether a mechanical part could be located at a particular position in the plane in order to facilitate manufacturing and maintenance processes [START_REF] Amundarain | Virtual reality for aircraft engines maintainability[END_REF] [65] [START_REF] Sagardia | Fast and accurate distance, penetration, and collision queries using point-sphere trees and distance fields[END_REF].

The main difficulty associated to haptic peripherals (those with force feedback) is that they need, to provide with a realistic sense of touch, a feedback response in the order of 500 Hz to 1 kHz. What this means in practice is that, very much like some 25 frames per second are needed in cinemas to provide the spectator with a continuous sensation of movement, in haptic environments nearly one thousand simulations per second must be carried out to provide the user with a continuous sense of touch. 1 kHz is roughly the free-hand gesture frequency [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF].

The reader will readily understand the difficulties associated with the simulation of non-linear solids (soft living tissues are frequently assumed to be hyperelastic, possibly with fiber reinforcement [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] [START_REF] Alastrue | Biomechanical modeling of refractive corneal surgery[END_REF]) under such astringent requirements. If, in addition, the possibility of contact between deformable solids is taken into account, the 1kHz constraint is even more difficult to fulfill.

Contact mechanics simulations under real-time restrictions have been tackled from a variety of approaches. One of the earliest and most popular is that of constructing Bounding Volume Hierarchies (BVH) [START_REF] Teschner | Collision detection for deformable objects[END_REF], consisting in associating each node in a tree with a subset of the primitives (polygons, NURBS, etc.) defining the boundary of the object. Other approaches for real time include the use of stochastic contact detection [START_REF] Uno | The sensitivity of presence to collision response[END_REF], based on the assumption that the perceived realism of contact detection depends more of the real-time response than in the accuracy itself of the simulation.

One of the most popular families of real-time contact detection algorithms is based on the use of distance fields [START_REF] Bridson | Simulation of clothing with folds and wrinkles[END_REF] [START_REF] Barbič | Six-dof haptic rendering of contact between geometrically complex reduced deformable models[END_REF]. Distance fields (level sets) constitute a very convenient way of representing very intricate geometries for contact detection, but has been traditionally considered as non-apt for real-time contact simulation between deformable solids, since the distance field must be updated according to the deformation of one of the solids [START_REF] Teschner | Collision detection for deformable objects[END_REF].

If we restrict ourselves to the problem of real-time simulation of hyperelastic solids, several approaches have been accomplished in the literature.

Besides the obvious choice of considering a purely elastic material, but which renders very poor results in terms of visual realism in the presence of large strains, the first approaches considered multi-resolution methods [START_REF] Kim | Virtual surgery simulation for medical training using multi-resolution organ models[END_REF] [START_REF] James | Multiresolution green's function methods for interactive simulation of large-scale elastostatic objects[END_REF].

Very popular in the last years are the methods based on the use of explicit finite element implementations [START_REF] Miller | Total lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF], possibly based on parallelization by using general-purpose Graphics Processing Units (GPUs) [START_REF] Grand | Real-time nonlinear finite element computations on {GPU} -application to neurosurgical simulation[END_REF] [59] [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF]. Due to the inherent complexity of the objective, techniques based upon model order reduction methods have also recently reached some popularity. Among them, we can cite [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF] [52] [48] [62] [61] [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF]. In general, all these last references are based on the use of Proper Orthogonal Decomposition (POD) techniques [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] [41] [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF], known also as Principal Component Analysis (PCA). These techniques employ an statistical treatment of the results obtained for complete, similar problems to the one at hand to construct a set of global, Ritz-like, basis functions that are optimal, with respecto to some measure, for the already solved, complete problems. They are then used for the problem under consideration in the hope that, if the modifications with respect to the original problems are small, they will be also a good choice for it. This approach presents some drawbacks. For instance, the choice of similar, complete problems, after which to obtain the set of global basis, is not an easy task [START_REF] Amsallem | An Interpolation Method for Adapting Reduced-Order Models and Application to Aeroelasticity[END_REF]. In addition, model reduction methods loose many of their advantages if we deal with non-linear problems, since they need for the reconstruction of the full tangent stiffness matrix in order to obtain a consistent linearization of the problem. Although several approaches have been developed to deal with this difficulty, no definitive response has been found [START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF] [14] [START_REF] Cueto | Real time simulation for computational surgery: a review[END_REF].

As an alternative to POD techniques, Proper Generalized Decomposition (PGD) techniques can be seen as a sort of a priori model order reduction method, in the sense that no complete problem must be solved in order to obtain the set of global basis functions. The origin of PGD dates back to the so-called radial loading approximation of the LArge Time INcrements (LATIN) method [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF]. It consisted, essentially, of an space-time separated representation of the solution that very much resembles that of POD, but obtained completely in an a priori fashion. Some years later, F. Chinesta employed a similar separated representation for the solution of high-dimensional problems arising from the kinetic theory of complex, non-Newtonian flows [3] [4]. Once both approaches have been identified as belonging to one single method for the model order reduction of PDEs, the field of application of PGD has grown exponentially. Just to cite some of the most recent survey papers, the interested reader can consult, for instance, [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] Based on PGD, the authors have developed a series of computational vademecum approaches [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF] for different physical problems such as thermal control of industrial furnaces [START_REF] Ch | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: Application to pultrusion[END_REF] [25], dynamic, data-driven application systems (DDDAS) [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF], shape optimization [START_REF] Ammar | Parametric solutions involving geometry: A step towards efficient shape optimization[END_REF] or computational surgery [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF],

to name a few. In essence, a computational vademecum (from the latin vade mecum, "goes with me") is a sort of reference guide that helps engineers (or physicians) by abridging known solutions to given problems. One of the earliest and most known examples is that of Bernoulli [12]. This concept has been translated by the authors to nowadays computational mechanics by computing, off-line and once for life, parametric (and thus multi-dimensional) solutions for a given problem. This general, multi-dimensional solution (a sort of computational response surface or meta-model, if preferred) is then used (evaluated) on-line at very high feedback rates.

The main difficulty with computational vademecums is that they give rise to high-dimensional problems. It is well-known that, for mesh-based methods such as finite elements or finite differences, high-dimensional problems suffer from the so-called curse of dimensionality [START_REF] Laughlin | The theory of everything[END_REF], i.e., an exponential increase on the number of degrees of freedom with the number of dimensions of the space. This difficulty is at the very heart of the choice of PGD for such high-dimensional problems, since it assumes the solution as being expressed as a finite sum of separable functions.

In this paper a novel approach to the real-time simulation of contact between non-linear solids is presented. It is based on an intensive use of computational vademecums obtained off-line by means of PGD techniques.

The developed technique can be seen as a particular instance of distance field methods, in which a general, high-dimensional distance field will be pre-computed off-line by taking the collision position as a parameter in the model. Thus, the most restrictive characteristic of distance field techniques will therefore be overcome, since there will be no need for distance field updating with the deformation of the solid.

In order to fully understand the developed methodology, a brief review of how to construct a computational vademecum by PGD techniques will be made in Section 2. The novel distance-field method will be presented in Section 3. In Section 4 examples of the performance of the proposed technique will be given. The paper is closed by some final remarks in Section 5.

A brief review of computational vademecums obtained by PGD techniques

As stated in the introduction, the main ingredient of the method here reported is the intensive use of computational vademecums [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF]. In essence, the idea behind this method is to pre-compute off-line every possible situation (with the limitations that will be later clarified) that could be faced during the on-line phase of the simulation.

Consider, as a starting point, the case of a vademecum in which we would like to store the displacement field u(x) of a non-linear (here, hyperelastic) solid Ω under the action of a force (assumed for the sake of simplicity as unitary and always acting in the vertical direction) at any point s of its boundary region Γ ⊂ Γ t ⊂ Γ = ∂Ω. This renders a problem defined in general in R 5 (u = u(x, s)), although if s is interpolated using nearestneighbour interpolation, it can be seen as a one-dimensional parameter (the node in which the load is acting), rendering a problem in R 4 .

In the general case of an arbitrary force acting in any direction, the problem would be stated so as to find the displacement field for any physical point

x ∈ R 3 , for a load acting at any point of its boundary (a two-dimensional manifold in R 3 ) and arbitrary values of its three components, thus t ∈ R 3 . This renders a problem defined in dimension 8.

Consider now the weak form of the equilibrium equations (balance of linear momentum). Again, for the sake of simplicity, we omit inertia terms.

For a detailed treatment of non-linear solid dynamics with the use of PGD and computational vademecums, the interested reader can consult [START_REF] Gonzalez | Real-time direct integration of reduced solid dynamics equations[END_REF]. Under these assumptions, the (doubly-) weak form of the problem, extended to the whole solid, Ω and the portion of its boundary which is accessible to load, Γ ⊂ Γ t , consists in finding the displacement u ∈ H 1 such that for all

u * ∈ H 1 0 : Γ Ω ∇ s u * : σdΩd Γ = Γ Γ t2 u * • tdΓd Γ (1)
where Γ = Γ u ∪ Γ t represents the boundary of the solid, divided into essential and natural regions, and where Γ t = Γ t1 ∪ Γ t2 , i.e., regions of homogeneous and non-homogeneous, respectively, natural boundary conditions. Here, t = -e k • δ(x -s), where δ represents the Dirac-delta function and e k the unit vector along the z-coordinate axis (we consider here, as mentioned before, and for the ease of exposition, a unit load directed towards the negative z axis of reference).

The Dirac-delta term is then regularized, and approximated by a truncated series of separable functions in the spirit of the PGD method, i.e.,

t j ≈ m i=1 f i j (x)g i j (s)
where m represents the order of truncation and f i j , g i j represent the j-th component of vectorial functions in space and boundary position, respectively. PGD techniques allow to efficiently construct the computational vademecum u(x, s) by constructing, in an iterative way, an approximation to the solution in the form of a finite sum of separable functions [START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF]. Let us assume that the method has converged to a solution, at iteration n of this procedure,

u n j (x, s) = n k=1 X k j (x) • Y k j (s), (2)
where the term u j refers to the j-th component of the displacement vector, j = 1, 2, 3 and functions X k (x) and Y k (s) represent the separated functions used to approximate the unknown field, obtained in previous iterations of the PGD algorithm. At this stage, the objetive of PGD is to provide the solution with an improvement given by the (n + 1)-th term of the approximation,

u n+1 j (x, s) = u n j (x, s) + R j (x) • S j (s), (3)
where R(x) and S(s) are the sought functions that improve the approximation. In an equivalent manner, admissible variations of this displacement field will be given by

u * j (x, s) = R * j (x) • S j (s) + R j (x) • S * j (s).
As can be noticed, even if the original problem is linear, PGD needs for the solution of a non-linear problem, i.e., to determine a product of functions, see Eq. (3). To this end, any of the possible linearization procedures could be employed. Traditionally, the authors have employed the simplest one, a fixed point, alternating directions algorithms that, although not being optimal, works in general very well. This strategy is briefly described hereafter.

Computation of S(s) assuming R(x) is known

In this case, following standard assumptions of variational calculus, we have

u * j (x, s) = R j (x) • S * j (s), (4)
or, equivalently, u * (x, s) = R • S * . The symbol "•" stands here for the so-called entry-wise, Hadamard or Schur multiplication for vectors. Once substituted into Eq. ( 1), gives

Γ Ω ∇ s (R • S * ) : C : ∇ s n k=1 X k • Y k + R • S dΩd Γ = Γ Γ t2 (R • S * ) • m k=1 f k • g k dΓd Γ,
or, equivalently (we omit obvious functional dependencies)

Γ Ω ∇ s (R • S * ) : C : ∇ s (R • S)dΩd Γ = Γ Γ t2 (R • S * ) • m k=1 f k • g k dΓd Γ - Γ Ω ∇ s (R • S * ) • R n dΩd Γ,
where R n represents:

R n = C : ∇ s u n .
Since the symmetric gradient operates on spatial variables only, we have:

Γ Ω (∇ s R • S * ) : C : (∇ s R • S)dΩd Γ = Γ Γ t2 (R • S * ) • m k=1 f k • g k dΓd Γ - Γ Ω (∇ s R • S * ) • R n dΩd Γ
where all the terms depending on x are known and hence all integrals over Ω and Γ t2 (the part of the natural boundary with non-homogeneous boundary conditions) can be computed to derive an equation to determine S(s).

Computation of R(x) assuming S(s) is known

Equivalently, in this case, we have

u * j (x, s) = R * j (x) • S j (s),
which, once substituted into Eq. ( 1), gives

Γ Ω ∇ s (R * • S) : C : ∇ s n k=1 X k • Y k + R • S dΩd Γ = Γ Γ t2 (R * • S) • m k=1 f k • g k dΓd Γ.
In this case all the terms depending on s (load position) can be integrated over Γ, leading to a generalized elastic problem to compute function R(x).

As mentioned before, even if fixed-point algorithms do not have guaranteed convergence properties, it is extremely unfrequent, in our experience, to find examples in which convergence is lost (see [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and references therein).

Other linearization strategies such as Newton-Raphson could also be equally employed [3] [4].

This development assumes implicitly small strain measures. For general, hyperelastic constitutive laws, large strain tensors (usually the Green-Lagrange tensor E) must be equally linearized. In the past, the authors have tackled this linearization problem within the PGD approach in two different ways. In [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF] an explicit approach was developed that renders, in general, very good results without stability problems. In [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] an approach was developed based on the combined use of PGD and Asymptotic Numerical Methods (ANM). In this last approach, the solution u is expanded in terms of a power series of an arc-length parameter, providing a sort of continuation method in which there is no need of updating tangent stiffness matrices. Other approaches such as the Discrete Empirical Interpolation Method [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] [14] could also be equally employed.

In general, following the aforementioned works, a strategy is followed in which enough terms in Eq. ( 2) so as to provide a norm of the residual below 10 -5 . Sometimes, given the strong limitations imposed by haptic response, further simplification of the PGD representation of the solution is needed.

In that cases, the number of terms are restricted to those giving an accuracy always above 10%, computed with respect to full FE solutions of the same problem.

In general, we refer the interested reader to the previous works in the field [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] [50] [START_REF] Gonzalez | Real-time direct integration of reduced solid dynamics equations[END_REF] for a thorough discussion on the relationship between accuracy and number of terms in the PGD representation of the solution.

A distance-field method based on the use of computational vedemecums

The method here developed assumes that a computational vademecum u(x, s) has already been computed for the solids under consideration. As mentioned before, this multi-dimensional solution provides a general solution for the displacement field in the solids under an arbitrary load acting on Γ.

The collision detection method assumes that one of the solids, say Ω 2 , is modeled as a pointshell [START_REF] Barbič | Six-dof haptic rendering of contact between geometrically complex reduced deformable models[END_REF], i.e., a set of boundary points (assumed for simplicity identical to the boundary nodes of the solid model, although this is not strictly necessary) equipped with normals to the surface. As such, the algorithm is not symmetric, since the choice of which body is chosen to be equipped with the pointshell affects the final result of the simulation, albeit slightly. The other solid, Ω 1 , is in turn equipped with a signed distance field, see Fig. 1. In our implementation, following closely [START_REF] Barbič | Six-dof haptic rendering of contact between geometrically complex reduced deformable models[END_REF], at every haptic cycle, contact penalty forces are determined by querying the points of Ω 2 against the distance field associated to Ω 1 . Traditionally, this approach has been considered not valid for haptic feedback requirements, if we deal with two deformable solids, see [START_REF] Teschner | Collision detection for deformable objects[END_REF]. This limitation is due to the need of updating the distance field along with the deformation of Ω 1 , which is not an easy task, even if the collision detection loop is not performed at every haptic cycle, as in [START_REF] Jeřábková | Stable cutting of deformable objects in virtual environments using xfem[END_REF]. In the approach here presented, there is no need for such an update, since a high-dimensional distance field is computed off-line that contains the distance fields for any deformed configuration of the solid.

Once collision has been detected, a force

F = -k c dn
is applied to both solids, provided that d is trilinearly interpolated from the distance field accompanying Ω 1 . k C is the contact penalty stiffness. The minus sign assumes that normals to Ω 2 point towards the interior of the solid, although this is completely arbitrary. Finally, n represents the normal to Ω 2 in the deformed configuration. An equal force is applied to Ω 1 , regardless of its geometry, since we have no information on its normals to the surface.

This very simple algorithm preserves continuity in the force computation, essential to perceive haptic response as realistic.

For very slender solids, large penetrations could eventually lead to decreasing force values, provided that Ω 2 crosses the medial axis of Ω 1 , see [START_REF] Barbič | Six-dof haptic rendering of contact between geometrically complex reduced deformable models[END_REF]. However, note that collision will be checked every millisecond. Except from very large solid velocities, this will likely never happen in fields such as computational surgery, for instance.

Remark 1. It is important to note that computational vademecums are de-signed for one single punctual load, or at least for a fixed number of loads. In contact problems the position of contact, and therefore the number of contacting nodes is a priori unknown. In the results presented in Section 4, where non-linear constitutive equations and strain measures are inherent, a predictor of the result is obtained by simply applying superposition, i.e., a linear combination of the contacting load predicted by the vademecums associated to each contacting node.

This very simple approach results to be enough for many applications such as computational surgery [START_REF] Cueto | Real time simulation for computational surgery: a review[END_REF], since visual realism is kept (no artificial gain of volume, as in purely elastic approaches, is obtained) and the human sense of touch is not able to detect the lack of accuracy generated in the resulting load transmitted by the haptic peripheral.

If more accuracy is needed, a correction can be added to this prediction by employing hyper reduction methods [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] [55] [56] [42]. In essence, the solution is projected onto the subspace spanned by the PGD modes of the solution (as in POD-based MOR) and the tangent stiffness matrix is evaluated only at a very limited number of elements of the mesh, typically some three times the number of modes considered. In our experiments, this approach has rendered excellent results, albeit not always necessary. This same rationale could be applied to the updating of the pre-computed distance field, although our experience indicates that the force feedback produced by the predictor obtained by linear superposition provides sufficient accuracy. I general, the correction is not perceivable by the human sense of touch.

The use of the correction step could seem compromising real time responses. However, it is important to notice that (i) only some elementary stiffness matrices are assembled (the global matrix is not singular because the global stiffness matrix appears pre-and post-multiplied by the matrix contain-ing the reduced basis) [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]; (ii) because the iteration starts from the prediction, that is in general close to the final solution, few correction iterations are needed; and (iii) only vector products are involved and all these can be computed simultaneously and in parallel by using GPUs, for example.

Even if the tests performed proved that surgical simulation applications do not need for the correction step, its consideration does not seem a great challenge, and it will be integrated in our future works to extend the applicability of this strategy to domains other than computational surgery.

Computation of the distance field

The main ingredient of this algorithm is the distance field associated to the solid Ω 1 . Any method could be in principle valid for the computation of this field. The key novelty of the approach here presented is that a highdimensional distance field

d = d(x, s)
for every load location s is computed off-line and stored in memory, very much in the spirit of computational vademecums. Here, x represents, by an abuse of notation, the coordinates of points belonging to the prism in which the distance field is computed, see Fig. 1. This prism, in general, does not need to be much bigger than the solid Ω 1 , see [START_REF] Teschner | Collision detection for deformable objects[END_REF], to avoid storing a big number of nodal distance values. However, in the implementation here presented, the distance field must be such that every possible deformation state of Ω 1 must lye within the prism. Once computed off-line and once for life, this distance field is evaluated for a particular load position s, that is, for any possible deformed configuration of the solid Ω 1 , thus giving a standard 3-d field for each load position, without the need of further updating.

Our implementation interpolates trilinearly the three physical coordinates and by nearest neighbors the s coordinate. This leads actually to a 4-d dis-tance field d = d(x, y, z, s), where s represents the node number to which the load is associated. Test codes used Matlab's griddedInterpolant algorithm [START_REF][END_REF] for a very fast evaluation of the distance field at a query point. Results

are given in Section 4 on the performance of this technique. The distance field at nodal position of the grid could be computed by employing fast marching algorithms [START_REF] Sethian | Fast marching methods[END_REF] or by simply employing the pdist Matlab function.

Remark 2. The brute force computation of the distance field is a feasible option if the number of dimensions of the problem (i.e., the number of parameters of the solution) is kept relatively small. If this is not so, the number of distance values to compute grows exponentially with the number of dimensions. The authors are currently exploring the possibility of solving the fast marching equations [START_REF] Sethian | Fast marching methods[END_REF] in parallel with the computation of the vademecum, Eqs. ( 2) and (3) in a PGD, separated way. In this manner, a vademecum could be obtained containing all the possible deformed configurations, each one attached to its distance field, both expressed in a separated form.

Numerical examples of performance

In what follows some examples have been selected to demonstrate the performance of the proposed method. Firstly, a Matlab [START_REF][END_REF] implementation of the algorithms was accomplished in order to check their validity. Once their performance had been assessed, a final implementation for the Geomagic Touch haptic device [START_REF] Geomagic | OpenHaptics Toolkit. 3D systems -Geomagic solutions[END_REF] was developed in order to test realism of the reaction forces, stability of the perceived touch, etc.

Contact between a beam and a rigid solid

To begin with, the simplest case of collision detection between a hyperelastic (Saint Venant Kirchhoff) cantilever beam and a solid is considered. A squared cross-section beam, with 40×40×400 mm, discretized with some 189 nodes, was considered. Young's modulus was assumed to be 209000 MPa, while Poisson coefficient was set to 0.3. The beam is actuated by a vertical load of 10 6 N, that can be applied at any point of the boundary of the beam.

If nearest-neighbor interpolation for the position of the load was assumed, a total of 170 different positions of this load could be encountered during the on-line phase of the simulation.

During the off-line phase of the calculation, a multidimensional solution is obtained for the displacement field of the beam, whose two first pairs of modes are shown in Fig. 2. Details of the PGD procedure for the solution of this parametric problem can be found at [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] or [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF]. A total of 81 pairs of functions X i and Y i (see Eq. ( 2)) were employed in this example, i.e., n = 81. Tests with only some 25 pairs also gave satisfactory results, with no perceptible loss of accuracy for the human eye. The off-line computation of this example takes roughly less than half an hour on a standard laptop.

Over this general solution u(x, s), a four-dimensional distance field d(x, s) is constructed. In this case, although the obvious choice would have been to

equip the rigid solid with the distance field, to demonstrate the capabilities of the method the inverse solution has been preferred. Therefore, the distance field was attached to the deformable solid, the beam. The mesh to construct this distance field contained 42 × 26 × 61 × 170 nodes (i.e., more than eleven millions of degrees of freedom) along x, y, z and s (load position, considered as a one-dimensional array of nodal locations) directions, respectively.

A rigid block is also considered with which the beam gets into contact, see Fig. running Matlab [START_REF][END_REF] the results indicated that the time to contact detection was 0.0006223 ± 1.02122E -07 seconds, well within the strong requirements of haptic contact feedback.

Contact between two hyperelastic beams

The fact that in the previous example a rigid solid was considered has no special implications on the complexity of the algorithm. If a deformabledeformable contact between two beams is now considered, the numbers remain roughly the same, the complexity of the algorithm depending essentially on the number of nodes of the pointshell and the number of degrees of freedom considered for the distance field. Both vademecums are essentially the same, so that the computer cost of the off-line computations remains the same of the previous example.

The algorithm here presented is so powerful that it can be implemented on an html page running javascript, for instance, and still provide visual perception of interactivity when the beam is touched with the mouse on the screen, see Fig. equipped with a Geomagic Touch [START_REF] Geomagic | OpenHaptics Toolkit. 3D systems -Geomagic solutions[END_REF] haptic device also gave excellent results, with no noticeable jumps in the perceived contact force, see Fig. 5.

Contact between two Stanford bunnies

In this example two Stanford bunnies were forced to contact each other in order to test the capabilities of the proposed technique. From the original STL file, a tetrahedral finite element model composed by 4285 nodes for each bunny was constructed, see Fig. 6. The total number of nodes in the pointshell (blue bunny in Fig. 6) was 3733. Bunnies were assumed to be modeled by a Kirchhoff-Saint Venant hyperelastic law with E = 2.0 • 10 7

MPa and ν = 0.3. The vademecum for each bunny -they are identical for both bunniescontained 143 modes. The cost of the off-line computation of these vademe-cums was more than three hours on a Mac Pro computer running Matlab [START_REF][END_REF] and equipped with a 6-core Intel Xeon ES CPU (although no parallelization was accomplished) at 3.5 GHz and with 16 GB RAM, 1867 MHz DDR3.

The red bunny is actuated by the user, provoking a rigid-solid displacement until contact occurs with the blue one. Three different positions were tested. At position number 1, no collision was detected, and the total time employed by the proposed algorithm was 0.0007148 seconds. At position number 2, 36 contacting nodes were detected, while the total time employed for their location was 0.001611 seconds. Finally, at position number three, 427 nodes got into contact, and in this last case 0.0012811 seconds were necessary to detect them.

The deformed configuration of the bunnies after contact is depicted in Fig. 7.

Again, for complex, detailed models of intricate geometry, the proposed method is able to render feedback at rates on the order of 1 kHz, proving the validity of the assumptions made so far.

Conclusions

In this paper a method for the real-time detection of collision between deformable solids has been presented. Although inspired by existing algorithms, the novelty of the proposed technique lies in the use of computational vademecums, i.e., pre-computed solutions that are exploited on-line at very high feedback rates. Here, the proposed method has been successfully tested on a Geomagic Touch device needing for 1 kHz feedback rates.

These extremely efficient vademecums are previously computed with the help of Proper Generalized Decomposition strategies, in order to alleviate the burden associated with the curse of dimensionality, i.e., the exponential The main advantage of the proposed method over existing algorithms is that it avoids the need of updating the distance field with the deformation of the reference solid. Since any deformation of the solid has been parametrized

previously by PGD and stored efficiently as a set of one-dimensional vectors, with the distance function associated to each deformed configuration, no need for on-line updating of the distance field results. This allows for a more detailed mesh for the computation of the distance field or, equivalently, for a more detailed geometric approximation of the pointshell in the second solid.

So far only collision detection has been accomplished within this algorithm, although be strongly believe that a generalization of the proposed technique is possible for frictional contact. This constitutes nowadays our main effort of research, whose results will be presented elsewhere.

  [START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF] [26] [16] [38] [29] [53].

Figure 1 :

 1 Figure 1: Sketch of the collision detection algorithm. Solid Ω 1 is equipped with a distance field, represented in the Figure, while solid Ω 2 is represented, for collision detection purposes, as a pointshell, i.e., a collection of boundary nodes and normals to the surface (not represented for simplicity).

Figure 2 :

 2 Figure 2: First two modes of the beam model, depending on the physical coordinates x (top) and load position s (bottom).

Figure 3 :

 3 Figure 3: Collision between a beam and a rigid block for one of the 170 possible load locations in the model. The beam (red), the block (blue) and a slice of the distance field are depicted.

4 .

 4 Contact iterations were limited to run in less than 100 ms, but the typical time to find collision was around 40 ms, still under the typical 25 frames per second of films.

Figure 4 :

 4 Figure 4: Collision between two beams running on an html web page with embedded javascript. The red stroke indicates the position of the load actuating the blue beam.

Figure 6 :

 6 Figure 6: Collision between two Stanford bunnies. A slice of the distance field is also shown.

Figure 7 :

 7 Figure 7: Deformed configuration of the bunnies after contact. Only the undeformed configuration of the blue bunny is represented, for simplicity, in white.

  [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF]. All the possible load cases were tested and the time necessary to detect contact stored. On a Macbook Pro laptop equipped with an Intel

	Z		Z
		Y		Y
	X	F	X	F
		0.4		0.3
		0.38		0.25
		0.36		0.2
		0.34		0.15
		0.32		0.1
		0.3		0.05
		0.28		0
		0.26	-0.05
		0.24	-0.1
		0.22	-0.15
		0.2	-0.2
		0.18	-0.25
		0.16	-0.3
		0.14	-0.35
		0.12	-0.4
		0.1	-0.45
		0.08	-0.5
		0.06
		0.04
		0.02
	(a)	(b)
	Z		Z
		Y		Y
	X	F	X	F
		-0.02	0.014
		-0.04	0.012
		-0.06	0.01
		-0.08	0.008
		-0.1	0.006
		-0.12	0.004
		-0.14	0.002
		-0.16	0
		-0.18	-0.002
		-0.2	-0.004
		-0.22	-0.006
		-0.24	-0.008
	(c)	(d)
	Core i7 processor running at 3 GHz (8Gb DDR3 RAM at 1600 MHz) and

$ This work has been partially supported by the Spanish Ministry of Science and Competitiveness, through grant number CICYT-DPI2011-27778-C02-01/02. Professor Chinesta is also supported by the Institut Universitaire de France. CIBER-BBN is an initiative funded by the VI National R&D&i Plan 2008-2011, Iniciativa Ingenio 2010, Consolider Program, CIBER Actions and financed by the Instituto de Salud Carlos III with assistance from the European Regional Development Fund.