1. Presentation of the problem
A common pitfall when solving problems with multiple physics is that the use of multiple tools often results in the creation of several meshes (see figure below). It is then a question of projecting the fields on a common grid in order to be able to carry out calculations at the same point in the space for the two calculated quantities.
In the following section, a mathematical tool is presented to project one field onto another, regular to no.
2. Projection methods
There are several algorithms of varying complexity to perform a field projection. Some variants are presented in this section.
2.1. Nearest neighbors
One technique for projecting a field onto another is to solve a minimization problem. The brief method is explained here given in code Aster in dimension two for the algorithm of the closest neighbors.
The aim is to project a cloud of points N bearing the field \varphi on a number of points M different of NLet us call \psi la projection de \varphi on M.
For each point k belonging M with coordinates (x_k,y_k), we're looking for its n nearest neighbors belonging to N Nas shown in following Figure.
The value of the field on the projected point \psi_k is computed with coefficients a_k,b_k,c_k such that:
\psi_k = a_k + b_k x_k + c_k y_k
It is therefore a question of finding, for all points k of cloud Mthe arrays [a_k], [b_k], [c_k] of coefficients that minimize the distance between the projected field \psi and the original field \varphiThus, for all point k of cloud M on cherche à minimiser les moindres carrés \varepsilon_k entre la fonction \psi_k and the original field \varphi_i on the n nearest neighbours
\varepsilon_k = \sum_{i=1}^{n} w(d_i) \big(\psi_k-\varphi_i \big)^2
Where the weighting function w(d_i) depends on distance d_i between point (x_i,y_i) of the initial cloud and point k of the projection grid (x_k,y_k) :
w(d_i) = e^{-(\frac{d_i}{d_{r}})^\beta}
Where
– d_i is the distance between (i-k) , i. e. d_i=\sqrt{(x_i-x_k)^2+(y_i-y_k)^2},
– d_r is the reference distance of the mesh. This is defined as the radius of the smallest circle centred at a point of M and containing three points of the cloud Nas shown in the figure below (this radius can be determined with a kdtree du package scipy par exemple). En première approximation, si la distance \Delta_x entre deux points du nuage est connue (cas d’un maillage régulier par exemple), on peut prendre d_r \sim \frac{\Delta x}{2}. Autrement, on pourra aussi utiliser la distance d_r \sim \frac{D^2}{N} = \frac{2D}{\sqrt{N}} where D est la distance entre les deux points les plus éloignés du nuage.
The choice of the parameter set (n,M,\beta,d_{r}) influe sur la performance globale de l’algorithme : selon la distance entre points du maillage il convient d’avoir un nombre de points de projection M while keeping in mind that the calculation time is proportional to M. What is more, chosing \beta>1implies that the influence of distant points is quickly amortized when their distance from the reference point is greater than d_{r}.
Les tests menés sur champs analytiques avec un densité de points homogène ont montré cependant que l’augmentation du nombre de voisins n’affecte pas la qualité de la projection. Pour de tels champs l’utilisation de n \in [4,8] est conseillée, avec pour coefficient de fonction de pondération \beta=1.5.
2.2 La méthode de Shepard
(élaboré à partir des travaux de T. Klinyam au sein d’AREP)
Cette méthode est une forme de pondération inverse à la distance qui est une méthode d’interpolation spatiale. L’effort principal ici va consister à construire une fonction d’interpolation Q(x,y) pour le champ considéré, au lieu d’utiliser directement la valeur du champ.
Sur les champs analytiques étudiés, cette méthode est \sim 2 \times plus performante en termes de précision, pour un temps d’exécution augmentant d’un ordre de grandeur. Si le gain peut paraître faible par rapport au temps de calcul, la robustesse de la méthode de Shepard pour la projection de champs réels a été éprouvée, notamment car elle induit moins de lissage sur ceux-ci.
Construction d’une fonction d’interpolation sur le champ initial N
La forme la plus simple de la méthode est appelée « Méthode de Shepard » (Shepard 1968) avec l’équation suivante :
f\left(x,y\right)=\sum_{k=1}^{N}{w_k\left(x,y\right) F_k}
Where w_k est une fonction de pondération généralement dépendante de la distance et F_k est une fonction qui permet de calculer la valeur au point k. Dans les cas les plus simples, comme la méthode des plus proches voisins, on définit F_k=\varphi_k, soit la valeur du champ initial au point k.
Dans la méthode de Shepard, on définit la pondération telle que :
w_k=\frac{\left[\frac{\left(R_q-d_k\right)_+}{R_q\cdot d_k}\right]^2}{\sum_{k=1}^{N}\left[\frac{\left(R_q-d_k\right)_+}{R_q\cdot d_k}\right]^2}
Where
– d_k est la distance entre le point (x,y) and point (x_k,y_k) , i. e. { d_k=\sqrt{\left(x-x_k\right)^2+\left(y-y_k\right)^2}}, dont \left(x,y\right) sont les coordonnées du point interpolé et (x_k,y_k) sont des coordonnées du nuage des points du champ de départ
– R_q est le rayon d’influence d’un cercle centré sur le point (x,y) (représenté sur la figure ci-dessous), En pratique on choisit N_q \sim 40 - 50. Le rayon d’influence est calculé de sorte que :
(R_q-d_k)_+= R_q - d_k \ \ \text{ pour }\ \ d_k < R_q
(R_q-d_k)_+=0 \text{ pour } d_k \ge R_q
Soit D la distance entre 2 points les plus éloignés du nuage N. On définit R_q via le paramètre de fonction nodale N_q such that
R_q=\frac{D}{2}\sqrt{\frac{N_q}{N}}
Pour ce qui est de la valeur de la fonction F_k, on construit ensuite une loi d’interpolation Q_k(x,y) définie de la manière suivante :
Q_k(x,y)=\varphi_k+a_{k2}(x-x_k)+a_{k3}(y-y_k)+a_{k4}(x-x_k)^2+a_{k5}(x-x_k)(y-y_k)\\+a_{k6}(y-y_k)^2
Where \varphi_k est la valeur sur le point (x_k,y_k) et les coefficients (a_{k2}, ...,\ a_{k6}) sont propres à chaque point (x_k,y_k). On les détermine par la minimisation des moindres carrés suivants, pour k \in [1,N^*], where N^* est le nombre de points dans le rayon d’influence R_q (afin de ne pas surcharger les notations on écrira indifféremment N^* and N dans la suite, en effet la fonction de pondération annule la contribution des autres points du nuage) :
\varepsilon_k = \sum_{ i=1\ et\ i\neq k}^{N^*}\left [ \frac{Q_k(x_i,y_i)-\varphi_i}{\rho_k(x_i,y_i)} \right ]^2
Soit, explicitement :
\varepsilon_k = \sum_{ i=1\ et\ i\neq k}^{N^*} \bigg [ \frac{(R_q-d_i)_+}{R_qd_i} \bigg( \varphi_k+a_{k2}(x_i-x_k)+a_{k3}(y_i-y_k)+a_{k4}(x_i-x_k)^2 \\ +a_{k5}(x_i-x_k)(y_i-y_k)+a_{k6}(y_i-y_k)^2 - \varphi_i \bigg) \bigg ]^2
With d_i la distance entre les points i and k : d_i= \sqrt{(x_k-x_i)^2+(y_k-y_i)^2}
L’équation d’interpolation devient alors :
f(x,y)=\sum_{k=1}^{N^*}{w_k(x,y) \cdot Q_k(x,y)}
On notera que l’effort calculatoire est situé dans l’étape d’identification des paramètres a_{2k},..., a_{6k} : plus le nombre N_q est grand, plus le nombre de points concerné est important, ce qui augmente la quantité de coefficients à déterminer pour la fonction de minimisation.
Calcul de la projection sur le champ M
On définit alors le rayon d’influence R_w, dont une illustration est donnée sur la Figure ci-après, tel que :
R_w=\frac{D}{2}\cdot \sqrt{\frac{N_w}{N}}
En pratique, un choix efficace consiste à prendre N_w=N_q/2.
La projection sur le point (x_m,y_m) of cloud M est calculée en faisant la somme des fonctions d’interpolations pour chaque point dans le rayon d’influence R_w (voir figure ci-après) :
\psi_m=\sum_{k=1}^{N}{w_k(x_m,y_m)}\cdot Q_k(x_m,y_m)
Il s’agit ensuite de calculer la projection \psi_m avec la fonction de pondération w_k définie de manière analogue à l’équation \eqref{eq:shep1_wk2} dont le rayon d’influence R_w est défini selon l’équation \eqref{eq:shep_rq} en utilisant le paramètre d’influence N_w.
Un exemple d’application de la projection selon la méthode de Shepard est donné ci-dessous.
Résultats sur le temps d’exécution des différents algorithmes de minimisation
On compare ici quelques unes des fonctions proposées la librairie optimize de scipy pour un problème de projection de champ analytique en 2D. De manière intéressante, la méthode L-BFGS-B est la plus performante. À partir de plus de 2000 points sur lesquels on calcule la fonction d’interpolation Q(x,y), on constate une augmentation linéaire du temps de calcul.
