initial commit

.gitattributes

@@ -0,0 +1,35 @@
+*.7z filter=lfs diff=lfs merge=lfs -text
+*.arrow filter=lfs diff=lfs merge=lfs -text
+*.bin filter=lfs diff=lfs merge=lfs -text
+*.bz2 filter=lfs diff=lfs merge=lfs -text
+*.ckpt filter=lfs diff=lfs merge=lfs -text
+*.ftz filter=lfs diff=lfs merge=lfs -text
+*.gz filter=lfs diff=lfs merge=lfs -text
+*.h5 filter=lfs diff=lfs merge=lfs -text
+*.joblib filter=lfs diff=lfs merge=lfs -text
+*.lfs.* filter=lfs diff=lfs merge=lfs -text
+*.mlmodel filter=lfs diff=lfs merge=lfs -text
+*.model filter=lfs diff=lfs merge=lfs -text
+*.msgpack filter=lfs diff=lfs merge=lfs -text
+*.npy filter=lfs diff=lfs merge=lfs -text
+*.npz filter=lfs diff=lfs merge=lfs -text
+*.onnx filter=lfs diff=lfs merge=lfs -text
+*.ot filter=lfs diff=lfs merge=lfs -text
+*.parquet filter=lfs diff=lfs merge=lfs -text
+*.pb filter=lfs diff=lfs merge=lfs -text
+*.pickle filter=lfs diff=lfs merge=lfs -text
+*.pkl filter=lfs diff=lfs merge=lfs -text
+*.pt filter=lfs diff=lfs merge=lfs -text
+*.pth filter=lfs diff=lfs merge=lfs -text
+*.rar filter=lfs diff=lfs merge=lfs -text
+*.safetensors filter=lfs diff=lfs merge=lfs -text
+saved_model/**/* filter=lfs diff=lfs merge=lfs -text
+*.tar.* filter=lfs diff=lfs merge=lfs -text
+*.tar filter=lfs diff=lfs merge=lfs -text
+*.tflite filter=lfs diff=lfs merge=lfs -text
+*.tgz filter=lfs diff=lfs merge=lfs -text
+*.wasm filter=lfs diff=lfs merge=lfs -text
+*.xz filter=lfs diff=lfs merge=lfs -text
+*.zip filter=lfs diff=lfs merge=lfs -text
+*.zst filter=lfs diff=lfs merge=lfs -text
+*tfevents* filter=lfs diff=lfs merge=lfs -text

Dockerfile

@@ -0,0 +1,34 @@
+# read the doc: https://huggingface.co/docs/hub/spaces-sdks-docker
+# you will also find guides on how best to write your Dockerfile
+
+FROM continuumio/anaconda3:main
+
+WORKDIR /code
+COPY ./environment.yml /code/environment.yml
+
+RUN apt-get update && apt-get install -y libgles2-mesa-dev
+
+# Create the environment using the environment.yml file
+RUN conda env create -f /code/environment.yml
+
+# Set up a new user named "user" with user ID 1000
+RUN useradd -m -u 1000 user
+# Switch to the "user" user
+USER user
+# Set home to the user's home directory
+ENV HOME=/home/user \
+    PYTHONPATH=$HOME/app \
+    PYTHONUNBUFFERED=1 \
+    GRADIO_ALLOW_FLAGGING=never \
+    GRADIO_NUM_PORTS=1 \
+    GRADIO_SERVER_NAME=0.0.0.0 \
+    GRADIO_THEME=huggingface \
+    SYSTEM=spaces
+
+# Set the working directory to the user's home directory
+WORKDIR $HOME/app
+
+# Copy the current directory contents into the container at $HOME/app setting the owner to the user
+COPY --chown=user . $HOME/app
+
+CMD ["./run.sh"]

README.md

@@ -0,0 +1,23 @@
+---
+title: MMGP demo
+emoji: 🚀
+colorFrom: purple
+colorTo: blue
+sdk: docker
+pinned: false
+app_port: 7860
+---
+
+This space provides a demo of a experiment from [the paper](https://arxiv.org/abs/2305.12871):
+
+```
+@misc{casenave2023mmgpmeshmorphinggaussian,
+      title={MMGP: a Mesh Morphing Gaussian Process-based machine learning method for regression of physical problems under non-parameterized geometrical variability}, 
+      author={Fabien Casenave and Brian Staber and Xavier Roynard},
+      year={2023},
+      eprint={2305.12871},
+      archivePrefix={arXiv},
+      primaryClass={cs.LG},
+      url={https://arxiv.org/abs/2305.12871}, 
+}
+```

app.py

@@ -0,0 +1,246 @@
+import gradio as gr
+import pickle
+from datasets import load_dataset
+from plaid.containers.sample import Sample
+
+
+import numpy as np
+import pyrender
+from trimesh import Trimesh
+import matplotlib as mpl
+import matplotlib.cm as cm
+
+from utils_inference import infer
+
+
+import os
+# switch to "osmesa" or "egl" before loading pyrender
+os.environ["PYOPENGL_PLATFORM"] = "egl"
+
+
+hf_dataset = load_dataset("PLAID-datasets/AirfRANS_remeshed", split="all_samples")
+
+file = open('training_data.pkl', 'rb')
+training_data = pickle.load(file)
+file.close()
+
+train_ids = hf_dataset.description['split']['ML4PhySim_Challenge_train']
+
+out_fields_names = hf_dataset.description['out_fields_names']
+in_scalars_names = hf_dataset.description['in_scalars_names']
+out_scalars_names = hf_dataset.description['out_scalars_names']
+
+
+nb_samples = len(hf_dataset)
+
+# <h2><b><a href='https://arxiv.org/abs/2305.12871' target='_blank'><b>MMGP</b> demo on the <a href='https://huggingface.co/datasets/PLAID-datasets/AirfRANS_remeshed' target='_blank'><b>AirfRANS_remeshed dataset</b></b></h2>
+
+# <a href='https://arxiv.org/abs/2305.12871' target='_blank'><b>MMGP paper</b>,
+
+_HEADER_ = '''
+<h2><b>MMGP demo on the <a href='https://huggingface.co/datasets/PLAID-datasets/AirfRANS_remeshed' target='_blank'><b>AirfRANS_remeshed dataset</b></b></h2>
+'''
+
+_HEADER_2 = '''
+The model is already trained. The morphing is the same as the one used in the [MMGP paper](https://arxiv.org/abs/2305.12871),
+but is much less involved than the one used in the winning entry of the [ML4PhySim competition](https://www.codabench.org/competitions/1534/).
+The training set has 103 samples and is the one used in this competition (some evaluations are out-of-distribution).
+
+The inference takes approx 5 seconds, and is done from scratch (no precomputation is used during the inference when evaluating samples).
+This means that the morphing and the finite element interpolations are re-done at each evaluation.
+
+After choosing a sample id, please change the field name in the dropdown menu to update the visualization.
+'''
+
+
+def round_num(num)->str:
+    return '%s' % float('%.3g' % num)
+
+
+def compute_inference(sample_id_str):
+    sample_id = int(sample_id_str)
+
+    sample_ = hf_dataset[sample_id]["sample"]
+    plaid_sample = Sample.model_validate(pickle.loads(sample_))
+
+    prediction = infer(hf_dataset, sample_id, training_data)
+    reference = {fieldn:plaid_sample.get_field(fieldn) for fieldn in out_fields_names}
+
+    nodes = plaid_sample.get_nodes()
+    if nodes.shape[1] == 2:
+        nodes__ = np.zeros((nodes.shape[0],nodes.shape[1]+1))
+        nodes__[:,:-1] = nodes
+        nodes = nodes__
+
+    triangles = plaid_sample.get_elements()['TRI_3']
+    trimesh = Trimesh(vertices = nodes, faces = triangles)
+
+    file = open('computed_inference.pkl', 'wb')
+    pickle.dump([trimesh, reference, prediction], file)
+    file.close()
+
+
+    str__ = f"Training sample {sample_id_str}"
+
+    if sample_id in train_ids:
+        str__ += " (in the training set)\n\n"
+    else:        
+        str__ += " (not in the training set)\n\n"
+
+    str__ += str(plaid_sample)+"\n"
+
+
+    if len(hf_dataset.description['in_scalars_names'])>0:        
+        str__ += "\nInput scalars:\n"
+        for sname in hf_dataset.description['in_scalars_names']:
+            str__ += f"- {sname}: {round_num(plaid_sample.get_scalar(sname))}\n"
+
+    str__ += f"\nNumber of nodes in the mesh: {nodes.shape[0]}"
+    
+    return str__
+
+
+def show_pred(fieldn):
+
+    file = open('computed_inference.pkl', 'rb')
+    data = pickle.load(file)
+    file.close()
+
+    trimesh, reference, prediction = data[0], data[1], data[2]
+
+    ref = reference[fieldn]
+    pred = prediction[fieldn]
+    
+    norm = mpl.colors.Normalize(vmin=np.min(ref), vmax=np.max(ref))
+    cmap = cm.seismic#cm.coolwarm
+    m = cm.ScalarMappable(norm=norm, cmap=cmap)
+    vertex_colors = m.to_rgba(pred)[:,:3]
+
+    trimesh.visual.vertex_colors = vertex_colors
+
+    mesh = pyrender.Mesh.from_trimesh(trimesh, smooth=False)
+
+    # compose scene
+    scene = pyrender.Scene(ambient_light=[.1, .1, .3], bg_color=[0, 0, 0])
+    camera = pyrender.PerspectiveCamera( yfov=np.pi / 3.0)
+    light = pyrender.DirectionalLight(color=[1,1,1], intensity=1000.)
+
+    scene.add(mesh, pose=  np.eye(4))
+    scene.add(light, pose=  np.eye(4))
+
+    scene.add(camera, pose=[[ 1,  0,  0,  1],
+                            [ 0,  1,  0,  0],
+                            [ 0,  0,  1,  6],
+                            [ 0,  0,  0,  1]])
+
+    # render scene
+    r = pyrender.OffscreenRenderer(1024, 1024)
+    color, _ = r.render(scene)
+
+    return color
+
+
+def show_ref(fieldn):
+
+    file = open('computed_inference.pkl', 'rb')
+    data = pickle.load(file)
+    file.close()
+
+    trimesh, reference, prediction = data[0], data[1], data[2]
+
+    ref = reference[fieldn]
+    
+    norm = mpl.colors.Normalize(vmin=np.min(ref), vmax=np.max(ref))
+    cmap = cm.seismic#cm.coolwarm
+    m = cm.ScalarMappable(norm=norm, cmap=cmap)
+    vertex_colors = m.to_rgba(ref)[:,:3]
+
+    trimesh.visual.vertex_colors = vertex_colors
+
+    mesh = pyrender.Mesh.from_trimesh(trimesh, smooth=False)
+
+    # compose scene
+    scene = pyrender.Scene(ambient_light=[.1, .1, .3], bg_color=[0, 0, 0])
+    camera = pyrender.PerspectiveCamera( yfov=np.pi / 3.0)
+    light = pyrender.DirectionalLight(color=[1,1,1], intensity=1000.)
+
+    scene.add(mesh, pose=  np.eye(4))
+    scene.add(light, pose=  np.eye(4))
+
+    scene.add(camera, pose=[[ 1,  0,  0,  1],
+                            [ 0,  1,  0,  0],
+                            [ 0,  0,  1,  6],
+                            [ 0,  0,  0,  1]])
+
+    # render scene
+    r = pyrender.OffscreenRenderer(1024, 1024)
+    color, _ = r.render(scene)
+
+    return color
+
+
+def show_err(fieldn):
+
+    file = open('computed_inference.pkl', 'rb')
+    data = pickle.load(file)
+    file.close()
+
+    trimesh, reference, prediction = data[0], data[1], data[2]
+
+    ref = reference[fieldn]
+    pred = prediction[fieldn]
+    
+    norm = mpl.colors.Normalize(vmin=np.min(ref), vmax=np.max(ref))
+    cmap = cm.seismic#cm.coolwarm
+    m = cm.ScalarMappable(norm=norm, cmap=cmap)
+    vertex_colors = m.to_rgba(pred-ref)[:,:3]
+
+    trimesh.visual.vertex_colors = vertex_colors
+
+    mesh = pyrender.Mesh.from_trimesh(trimesh, smooth=False)
+
+    # compose scene
+    scene = pyrender.Scene(ambient_light=[.1, .1, .3], bg_color=[0, 0, 0])
+    camera = pyrender.PerspectiveCamera( yfov=np.pi / 3.0)
+    light = pyrender.DirectionalLight(color=[1,1,1], intensity=1000.)
+
+    scene.add(mesh, pose=  np.eye(4))
+    scene.add(light, pose=  np.eye(4))
+
+    scene.add(camera, pose=[[ 1,  0,  0,  1],
+                            [ 0,  1,  0,  0],
+                            [ 0,  0,  1,  6],
+                            [ 0,  0,  0,  1]])
+
+    # render scene
+    r = pyrender.OffscreenRenderer(1024, 1024)
+    color, _ = r.render(scene)
+
+    return color
+
+
+if __name__ == "__main__":
+
+    with gr.Blocks() as demo:
+
+        # trimesh, reference, prediction = compute_inference(0)
+        gr.Markdown(_HEADER_)
+        gr.Markdown(_HEADER_2)
+        with gr.Row(variant="panel"):
+            with gr.Column():
+                d1 = gr.Slider(0, nb_samples-1, value=0, label="Training sample id", info="Choose between 0 and "+str(nb_samples-1))
+                # output1 = gr.Text(label="Inference status")
+                output4 = gr.Text(label="Information on sample")
+                output5 = gr.Image(label="Error")      
+            with gr.Column():
+                d2 = gr.Dropdown(out_fields_names, value=out_fields_names[0], label="Field name")        
+                output2 = gr.Image(label="Reference")      
+                output3 = gr.Image(label="MMGP prediction")
+
+        # d1.input(compute_inference, [d1], [output1, output4])
+        d1.input(compute_inference, [d1], [output4])
+        d2.input(show_ref, [d2], [output2]) 
+        d2.input(show_pred, [d2], [output3]) 
+        d2.input(show_err, [d2], [output5]) 
+    
+    demo.launch()

environment.yml

@@ -0,0 +1,15 @@
+name: demo
+channels:
+  - conda-forge
+  - nodefaults
+
+dependencies:
+  - python=3.11
+  - muscat
+  - plaid
+  - datasets
+  - pyrender
+  - matplotlib
+  - gradio
+  - pytorch
+  # - pyopengl

gp_torch.py

@@ -0,0 +1,93 @@
+import numpy as np
+import torch
+import torch.nn as nn
+
+class GaussianProcessRegressor(nn.Module):
+    def __init__(
+        self,
+        length_scale=1.0,
+        noise_scale=1.0,
+        amplitude_scale=1.0,
+    ):
+        super().__init__()
+        if isinstance(length_scale, float):
+            length_scale = np.array([length_scale])
+        elif isinstance(length_scale, np.ndarray):
+            assert length_scale.ndim == 1
+        else:
+            raise TypeError()
+
+        self.register_parameter(
+            "length_scale_",
+            param=nn.Parameter(torch.Tensor(np.log(length_scale)), requires_grad=True),
+        )
+        self.register_parameter(
+            "noise_scale_",
+            param=nn.Parameter(torch.tensor(np.log(noise_scale)), requires_grad=True),
+        )
+        self.register_parameter(
+            "amplitude_scale_",
+            param=nn.Parameter(
+                torch.tensor(np.log(amplitude_scale)), requires_grad=True
+            ),
+        )
+
+        self.nll = None
+
+    def forward(self, x):
+        alpha = self.alpha
+        k = self.Kxy(self.X, x)
+        mu = k.T.mm(alpha)
+        return mu
+
+    def log_marginal_likelihood(self, X, y):
+        D = X.shape[1]
+        K = self.Kxx(X)
+        L = torch.linalg.cholesky(K)
+        alpha = torch.linalg.solve(L.T, torch.linalg.solve(L, y))
+        marginal_likelihood = (
+            -0.5 * y.T.mm(alpha)
+            - torch.log(torch.diag(L)).sum()
+            - D * 0.5 * np.log(2 * np.pi)
+        )
+        self.L = L
+        self.alpha = alpha
+        self.K = K
+        return marginal_likelihood
+
+    def Kxx(self, X):
+        param = self.length_scale_.exp().sqrt()
+        sqdist = torch.cdist(X / param[None], X / param[None]) ** 2
+
+        res = self.amplitude_scale_.exp() * torch.exp(-0.5 * sqdist) + self.noise_scale_.exp() * torch.eye(len(X)).type_as(X)
+
+        return res
+
+    def Kxy(self, X, Z):
+        param = self.length_scale_.exp().sqrt()
+        sqdist = torch.cdist(X / param[None], Z / param[None]) ** 2
+        res = self.amplitude_scale_.exp() * torch.exp(-0.5 * sqdist)
+
+        return res
+
+    def fit(self, X, y, opt, num_steps):
+        assert X.shape[1] == len(self.length_scale_)
+        self.y = y
+        self.X = X
+
+        scheduler = torch.optim.lr_scheduler.ExponentialLR(opt, gamma=0.9)
+
+        self.train()
+        nll_hist = []
+        for it in range(num_steps):
+            opt.zero_grad()
+            try:
+                nll = -self.log_marginal_likelihood(self.X, self.y).sum()
+            except torch.linalg.LinAlgError:
+                break
+            nll.backward()
+            opt.step()
+            if it%10==0 and it<1000:
+                scheduler.step()
+            nll_hist.append(nll.item())
+        return nll_hist

old_gradio_config/default_docker_config.txt

@@ -0,0 +1,58 @@
+# FROM python:3.10
+
+# RUN useradd -m -u 1000 user 
+
+# RUN apt-get update && apt-get install -y fakeroot &&     mv /usr/bin/apt-get /usr/bin/.apt-get &&     echo '#!/usr/bin/env sh\nfakeroot /usr/bin/.apt-get $@' > /usr/bin/apt-get &&     chmod +x /usr/bin/apt-get && 	rm -rf /var/lib/apt/lists/* && 	useradd -m -u 1000 user
+
+# COPY --chown=1000:1000 --from=root / /
+
+# WORKDIR /home/user/app
+
+# RUN apt-get update && apt-get install -y 	git 	git-lfs 	ffmpeg 	libsm6 	libxext6 	cmake 	rsync 	libgl1-mesa-glx  libgles2-mesa-dev	&& rm -rf /var/lib/apt/lists/* 	&& git lfs install
+
+# RUN pip install --no-cache-dir pip==22.3.1 && 	pip install --no-cache-dir 	datasets 	"huggingface-hub>=0.19" "hf-transfer>=0.1.4" "protobuf<4" "click<8.1" "pydantic~=1.0"
+
+# RUN --mount=target=/tmp/pre-requirements.txt,source=pre-requirements.txt 	pip install --no-cache-dir -r /tmp/pre-requirements.txt
+
+# RUN --mount=target=/tmp/requirements.txt,source=requirements.txt     pip install --no-cache-dir -r /tmp/requirements.txt
+
+# RUN export MUSCAT_USE_EIGENCYEIGEN=1
+
+# RUN wget https://boostorg.jfrog.io/artifactory/main/release/1.82.0/source/boost_1_82_0.zip
+
+# RUN unzip boost_1_82_0.zip
+
+# RUN MUSCAT_USE_EIGENCYEIGEN=1
+
+# RUN MUSCAT_EXTRA_INCLUDE_DIRS=%cd%\boost_1_82_0
+
+# RUN pip install Muscat@https://gitlab.com/drti/muscat/-/archive/2.2.1/muscat-2.0.2.tar.bz2
+
+# RUN pip install --no-cache-dir 	gradio[oauth]==4.37.1 	"uvicorn>=0.14.0" 	spaces
+
+# COPY --link --chown=1000 ./ /home/user/app
+
+# RUN pip freeze > /tmp/freeze.txt
+
+# COPY --from=pipfreeze --link --chown=1000 /tmp/freeze.txt /tmp/freeze.txt
+
+
+
+
+# # Set up a new user named "user" with user ID 1000
+# RUN useradd -m -u 1000 user
+
+# # Switch to the "user" user
+# USER user
+
+# # Set home to the user's home directory
+# ENV HOME=/home/user \
+#     PATH=/home/user/.local/bin:$PATH
+
+# # Set the working directory to the user's home directory
+# WORKDIR $HOME/app
+
+# # Copy the current directory contents into the container at $HOME/app setting the owner to the user
+# COPY --chown=user . $HOME/app
+
+# CMD ["python", "app.py"]

old_gradio_config/packages.txt

@@ -0,0 +1 @@
+libgles2-mesa-dev

old_gradio_config/pre-requirements.txt

@@ -0,0 +1,9 @@
+eigency
+mkl
+numpy
+sympy
+mkl-include
+cython
+wheel
+tbb-devel
+dill

old_gradio_config/requirements.txt

@@ -0,0 +1,9 @@
+git+https://gitlab.com/drti/plaid@version_0.0.10_python310
+numpy==1.26.4
+tqdm
+pyyaml
+pycgns
+scikit-learn
+datasets
+pyrender
+matplotlib

run.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+CONDA_ENV=$(head -1 /code/environment.yml | cut -d" " -f2)
+eval "$(conda shell.bash hook)"
+conda activate $CONDA_ENV
+python app.py

training_data.pkl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5a70d30164ee7489c55edb72d37a48c78a62a9d846b2fd2bda55574984537244
+size 22328714

utils_inference.py

@@ -0,0 +1,380 @@
+import Muscat.Containers.ElementsDescription as ED
+
+from Muscat.FE.FETools import PrepareFEComputation
+from Muscat.FE.Fields.FEField import FEField
+import Muscat.Containers.MeshModificationTools as MMT
+
+
+from Muscat.Containers import MeshGraphTools as MGT
+from Muscat.Containers.Filters import FilterObjects as FO
+from Muscat.Containers.Filters import FilterOperators as FOp
+import Muscat.Containers.MeshInspectionTools as MIT
+
+
+out_fields_names = ['Ux', 'Uy', 'p', 'nut']
+in_scalars_names = ['angle_of_attack', 'inlet_velocity']
+
+from Muscat.Containers.NativeTransfer import NativeTransfer
+import copy
+import bisect
+import numpy as np
+from scipy import spatial
+import pickle
+from Muscat.Bridges.CGNSBridge import CGNSToMesh
+
+
+import torch
+device = torch.device("cpu")
+
+
+num_steps = [2_000]
+
+from plaid.containers.dataset import Sample
+
+
+
+
+
+def pretreat_sample(mesh):
+    ###################
+    # Compute the skin of the mesh (containing the external boundary and the airfoil boundary)
+    ###################
+    MMT.ComputeSkin(mesh, md = 2, inPlace = True)
+
+    ###################
+    # Extract ids of the bar elements corresponding to the airfoil
+    ###################
+    ff1 = FO.ElementFilter(zone = lambda p: (-p[:,0]-1.99))
+    ff2 = FO.ElementFilter(zone = lambda p: (p[:,0]-3.99))
+    ff3 = FO.ElementFilter(zone = lambda p: (-p[:,1]-1.49))
+    ff4 = FO.ElementFilter(zone = lambda p: (p[:,1]-1.49))
+    efAirfoil = FOp.IntersectionFilter(filters=[ff1, ff2, ff3, ff4])
+    airfoil_ids = efAirfoil.GetIdsToTreat(mesh, "bar2")
+
+    ###################
+    # Preparations
+    ###################
+    ext_bound = np.setdiff1d(mesh.elements["bar2"].GetTag("Skin").GetIds(), airfoil_ids)
+    mesh.elements["bar2"].GetTag("External_boundary").SetIds(ext_bound)
+    mesh.elements["bar2"].GetTag("Airfoil").SetIds(airfoil_ids)
+    nfExtBoundary = FO.NodeFilter(eTag = "External_boundary")
+    nodeIndexExtBoundary  = nfExtBoundary.GetNodesIndices(mesh)
+    mesh.GetNodalTag("External_boundary").AddToTag(nodeIndexExtBoundary)
+    nfAirfoil = FO.NodeFilter(eTag = "Airfoil")
+    nodeIndexAirfoil  = nfAirfoil.GetNodesIndices(mesh)
+    mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil)
+
+
+    ###################
+    # Add node tag for the intrado and extrado
+    ###################
+    nfAirfoil = FO.NodeFilter(eTag = "Airfoil")
+    nodeIndexAirfoil = nfAirfoil.GetNodesIndices(mesh)
+    mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil)
+
+    airfoil = ExtractAirfoil(mesh)
+    indices_extrado = airfoil[0][0]
+    indices_intrado = airfoil[0][1]
+    mesh.GetNodalTag("Extrado").AddToTag(indices_extrado)
+    mesh.GetNodalTag("Intrado").AddToTag(indices_intrado)
+
+    efExtrado = FO.ElementFilter(nTag = "Extrado")
+    efIntrado = FO.ElementFilter(nTag = "Intrado")
+    mesh.elements["bar2"].GetTag("Extrado").SetIds(efExtrado.GetIdsToTreat(mesh, "bar2"))
+    mesh.elements["bar2"].GetTag("Intrado").SetIds(efIntrado.GetIdsToTreat(mesh, "bar2"))
+
+
+def ExtractPathFromMeshOfBars(mesh, startingClosestToPoint, trigo_dir = True):
+
+    nodeGraph0Airfoild = MGT.ComputeNodeToNodeGraph(mesh, dimensionality=1)
+    nodeGraphAirfoild = [list(nodeGraph0Airfoild[i].keys()) for i in range(nodeGraph0Airfoild.number_of_nodes())]
+
+    tree = spatial.KDTree(mesh.nodes)
+    _, indicesTrailEdge = tree.query([startingClosestToPoint], k=1)
+
+    p1init = indicesTrailEdge[0]
+
+    temp1=mesh.nodes[nodeGraphAirfoild[p1init][0]][1]
+    temp2=mesh.nodes[nodeGraphAirfoild[p1init][1]][1]
+
+    if trigo_dir:
+        condition = temp1 > temp2
+    else:
+        condition = temp1 < temp2
+
+    if condition:
+        p2 = nodeGraphAirfoild[p1init][0]
+    else:
+        p2 = nodeGraphAirfoild[p1init][1]
+
+    p1 = p1init
+    path = [p1, p2]
+    while p2 != p1init:
+        p2save = p2
+        tempArray = np.asarray(nodeGraphAirfoild[p2])
+        p2 = tempArray[tempArray!=p1][0]
+        p1 = p2save
+        path.append(p2)
+
+    return path
+
+
+def ExtractAirfoil(mesh):
+
+    efAirfoil = FO.ElementFilter(elementType=ED.Bar_2, eTag=["Airfoil"])
+    airfoilMesh = MIT.ExtractElementsByElementFilter(mesh, efAirfoil)
+
+    path = ExtractPathFromMeshOfBars(airfoilMesh, np.array([1.,0.]))
+
+    tree = spatial.KDTree(airfoilMesh.nodes[path])
+    _, indicesLeadEdge = tree.query([[0.,0.]], k=1)
+
+    indices_extrado = path[:indicesLeadEdge[0]+1]
+    indices_intrado = path[indicesLeadEdge[0]:]
+
+    indices_airfoil = [indices_extrado, indices_intrado]
+
+    nodes_extrado = mesh.nodes[indices_extrado]
+    nodes_intrado = mesh.nodes[indices_intrado]
+
+    nodes_airfoil = [nodes_extrado, nodes_intrado]
+
+    return indices_airfoil, nodes_airfoil
+
+
+def computeAirfoilCurvAbscissa(airfoil):
+
+    indices_airfoil = airfoil[0]
+    nodes_airfoil = airfoil[1]
+
+    curv_abscissa = []
+    for i in range(2):
+        local_curv_abscissa = np.zeros(len(indices_airfoil[i]))
+        for j in range(1,len(local_curv_abscissa)):
+            local_curv_abscissa[j] = local_curv_abscissa[j-1] + np.linalg.norm(nodes_airfoil[i][j]-nodes_airfoil[i][j-1])
+        local_curv_abscissa /= local_curv_abscissa[-1]
+        curv_abscissa.append(local_curv_abscissa)
+
+    return curv_abscissa
+
+
+def MapAirfoil(airfoil_ref, curv_abscissa_ref, curv_abscissa):
+
+    nodes_airfoil_ref = airfoil_ref[1]
+    dim_nodes = nodes_airfoil_ref[0][0].shape[0]
+
+    mapped_airfoil = []
+    for i in range(2):
+        local_mapped_airfoil = np.zeros((len(curv_abscissa[i])-1, dim_nodes))
+        for j in range(len(curv_abscissa[i])-1):
+            index = max(bisect.bisect_right(curv_abscissa_ref[i], curv_abscissa[i][j]) - 1, 0)
+
+            a = nodes_airfoil_ref[i][index]
+            b = nodes_airfoil_ref[i][index+1]
+            dl = curv_abscissa[i][j] - curv_abscissa_ref[i][index]
+            dir = (b-a)/np.linalg.norm(b-a)
+            local_mapped_airfoil[j] = a + dl * dir
+        mapped_airfoil.append(local_mapped_airfoil)
+
+    return mapped_airfoil
+
+
+
+
+def GetFieldTransferOpCppStep1(inputField, nbThreads = None):
+
+    method="Interp/Clamp"
+
+    nt = NativeTransfer()
+
+    if nbThreads is not None:
+        nt.SetMaxConcurrency(nbThreads)
+
+    nt.SetTransferMethod(method)
+    defaultOptions = {"usePointSearch": True,
+                    "useElementSearch": False,
+                    "useElementSearchFast": False,
+                    "useEdgeSearch": True,
+                    }
+
+    options = {}
+
+    defaultOptions.update(options)
+
+    dispatch = {"usePointSearch": nt.SetUsePointSearch,
+                "useElementSearch": nt.SetUseElementSearch,
+                "useElementSearchFast": nt.SetUseElementSearchFast,
+                "useEdgeSearch": nt.SetUseEdgeSearch,
+                "DifferentialOperator": nt.SetDifferentialOperator,
+                }
+
+    for k, v in defaultOptions.items():
+        if k in dispatch.keys():
+            dispatch[k](v)
+        else:
+            raise RuntimeError(f"Option {k} not valid")
+
+    nt.SetSourceFEField(inputField, None)
+
+    return nt
+
+
+
+def GetFieldTransferOpCppStep2(nt, targetPoints):
+
+    nt.SetTargetPoints(targetPoints)
+
+    nt.Compute()
+    op = nt.GetOperator()
+    status = nt.GetStatus()
+    return op, status
+
+
+
+
+def morph_sample(mesh, airfoil_0, curv_abscissa_0):
+
+    airfoil_1 = ExtractAirfoil(mesh)
+    curv_abscissa_1 = computeAirfoilCurvAbscissa(airfoil_1)
+
+    mapped_airfoil = MapAirfoil(airfoil_0, curv_abscissa_0, curv_abscissa_1)
+
+    ##############################################################
+    # Compute global target displacement and masks for RBF field morphing
+    ##############################################################
+    indices_extrado_to_morph_1 = airfoil_1[0][0][:-1]
+    indices_intrado_to_morph_1 = airfoil_1[0][1][:-1]
+    other_boundary_ids_1 = mesh.GetNodalTag("External_boundary").GetIds()
+
+    l1 = len(indices_extrado_to_morph_1)
+    l2 = len(indices_intrado_to_morph_1)
+    l3 = len(other_boundary_ids_1)
+
+    targetDisplacement     = np.zeros((l1 + l2 + l3, 2))
+    targetDisplacementMask = np.zeros((l1 + l2 + l3), dtype = int)
+
+    targetDisplacement[:l1]                        = mapped_airfoil[0] - mesh.nodes[indices_extrado_to_morph_1]
+    targetDisplacement[l1:l1+l2]                   = mapped_airfoil[1] - mesh.nodes[indices_intrado_to_morph_1]
+    targetDisplacement[l1+l2:l1+l2+l3]             = np.zeros((l3,2))
+
+    targetDisplacementMask[:l1]                        = indices_extrado_to_morph_1
+    targetDisplacementMask[l1:l1+l2]                   = indices_intrado_to_morph_1
+    targetDisplacementMask[l1+l2:l1+l2+l3]             = other_boundary_ids_1
+
+    ##############################################################
+    # Compute the morphing
+    ##############################################################
+
+    # RBF morphing
+
+    mesh_nodes = mesh.nodes.copy()
+    morphed_nodes = MMT.Morphing(mesh, targetDisplacement, targetDisplacementMask, radius=None)
+    mesh.nodes = morphed_nodes
+
+    mesh.nodeFields['X'] = mesh_nodes[:,0]
+    mesh.nodeFields['Y'] = mesh_nodes[:,1]
+
+    return mesh
+
+
+def project_sample(morphed_mesh, morphed_mesh_0):
+    projected_mesh = copy.deepcopy(morphed_mesh_0)
+
+    space_, numberings_, _, _ = PrepareFEComputation(morphed_mesh)
+    inputFEField = FEField(name="dummy", mesh=morphed_mesh, space=space_, numbering=numberings_[0])
+
+    nt = GetFieldTransferOpCppStep1(inputFEField, 1)
+    FE_interpolation_op, _ = GetFieldTransferOpCppStep2(nt, morphed_mesh_0.nodes)
+
+    for pfn in out_fields_names + ['X', 'Y']:
+        projected_mesh.nodeFields[pfn] = FE_interpolation_op.dot(morphed_mesh.nodeFields[pfn])
+
+    return projected_mesh
+
+
+
+
+def pretreat_morph_and_project_mesh(i_sample, dataset, indices, morphed_mesh_0, airfoil_0, curv_abscissa_0):
+
+
+    ###################
+    ## 1) Pretreat data
+    ###################
+    sample = Sample.model_validate(pickle.loads(dataset[int(indices[i_sample])]["sample"]))
+    mesh = CGNSToMesh(sample.get_mesh())
+    pretreat_sample(mesh)
+
+    ###################
+    ## 2) Morph data
+    ###################
+    morphed_mesh = morph_sample(mesh, airfoil_0, curv_abscissa_0)
+
+    # print("morphed_mesh =", morphed_mesh)
+
+    # from Muscat.IO import XdmfWriter as XW
+    # XW.WriteMeshToXdmf("morphed_mesh.xdmf", morphed_mesh)
+
+    ###################
+    ## 3) Project data
+    ###################
+    projected_mesh = project_sample(morphed_mesh, morphed_mesh_0)
+
+    # print("projected_mesh =", projected_mesh)
+
+    # from Muscat.IO import XdmfWriter as XW
+    # XW.WriteMeshToXdmf("projected_mesh.xdmf", projected_mesh)
+
+    in_scalars = [sample.get_scalar(sn) for sn in in_scalars_names]
+
+    return [projected_mesh, in_scalars, morphed_mesh.nodes]
+
+
+def infer(dataset, indice, training_data):
+
+    common_mesh, correlationOperator2c, airfoil_0, curv_abscissa_0, scalar_scalers, scalerX, pca_clouds, pca_fields, y_scalers, kmodels, X_train, y_train = training_data
+
+    projected_mesh, in_scalars, morphed_mesh_nodes = pretreat_morph_and_project_mesh(indice, dataset = dataset, indices = np.arange(len(dataset)), morphed_mesh_0 = common_mesh, airfoil_0 = airfoil_0, curv_abscissa_0 = curv_abscissa_0)
+
+    space_, numberings_, _, _ = PrepareFEComputation(common_mesh)
+    inputFEField_0 = FEField(name="dummy", mesh=common_mesh, space=space_, numbering=numberings_[0])
+    nt_0 = GetFieldTransferOpCppStep1(inputFEField_0, nbThreads = 2)
+
+    iffo = GetFieldTransferOpCppStep2(nt_0, morphed_mesh_nodes)[0]
+
+    clouds = np.stack([projected_mesh.nodeFields["X"], projected_mesh.nodeFields["Y"]], axis=1)
+
+    scalars = np.array(in_scalars)
+
+    clouds = clouds.reshape(-1, 1)
+
+    X_pca = np.dot(pca_clouds, correlationOperator2c.dot(clouds)).T
+    X_scalars = scalar_scalers.transform(scalars.reshape(1, -1))
+    unscaled_X = np.concatenate([X_pca, X_scalars], axis=-1)
+    X = scalerX.transform(unscaled_X)
+
+
+    predictions = {}
+    y = []
+    for i, fn in enumerate(out_fields_names):
+
+        X_ = torch.tensor(X, dtype=torch.float64).to(device)
+
+        output_dim = pca_fields[0].shape[0]
+
+        n_samples = 1
+        y_pred_i = np.empty((n_samples, output_dim, len(num_steps)))
+
+        for j in range(output_dim):
+            for k in range(len(num_steps)):
+                y_pred_i[:,j,k] = kmodels[i][j][k](X_).detach().cpu().numpy().squeeze()
+
+        y.append(y_pred_i)
+        y_pred_i = np.mean(y_pred_i, axis = 2)
+
+        y_pred_i = y_scalers[i].inverse_transform(y_pred_i)
+
+        y_pred_common_i = np.dot(y_pred_i, pca_fields[i]).flatten()
+
+        predictions[fn] = iffo.dot(y_pred_common_i)
+
+    return predictions