7  FEniCS-Link — Optimisation FEM par DAG

7.1 🎯 Objectif du TP

L’objectif de ce travail pratique dépasse la simple résolution d’un problème de mécanique.

Il s’agit ici de :

  • découvrir un environnement de calcul basé sur un DAG (Directed Acyclic Graph)
  • construire un pipeline scientifique complet : simulation → analyse → optimisation
  • automatiser une étude de convergence de maillage à l’aide d’une approche HPO (Hyperparameter Optimization)

7.2 🧭 Contexte pédagogique

Dans une approche classique en éléments finis : - on choisit un maillage - on exécute une simulation - on ajuste les paramètres manuellement (approche essai/erreur)

Cette démarche est : - coûteuse - peu reproductible - difficilement automatisable

7.3 🚀 Approche proposée

Dans ce TP, on adopte une approche moderne basée sur un pipeline de calcul structuré :

💡 Un workflow automatisé piloté par un DAG

Ce pipeline est exécuté via Link (extension JupyterLab), qui permet de :

  • structurer le notebook sous forme de graphe de calcul
  • gérer automatiquement les dépendances entre étapes
  • éviter les recalculs inutiles (caching)
  • lancer des optimisations de manière automatisée

7.4 🧪 Problème étudié

On considère une poutre encastrée soumise à une charge, utilisée comme cas test.

Ce problème permet d’illustrer :

  • la construction d’un modèle éléments finis
  • l’analyse des résultats mécaniques
  • la quantification de l’erreur
  • l’optimisation automatique du maillage

7.5 🧩 Vue globale du pipeline

Le pipeline mis en place suit les étapes suivantes :

  1. Chargement des paramètres du problème
  2. Génération du maillage (nx, ny, nz)
  3. Simulation FEM
  4. Post-traitement (déplacement, contraintes)
  5. Calcul de l’erreur
  6. Optimisation des paramètres de maillage”

7.6 ♻️ Interprétation en DAG

Dans Link, chaque étape devient un nœud du graphe :

Inputs → Mesh → FEM → Post → Erreur → HPO
                   ↑_____________________|

👉 La boucle introduite par le HPO permet d’automatiser la convergence.

7.7 🕸️ Illustration du pipeline

FEniCS-Link

Figure : DAG pour FEniCS

7.8 ⚡ Philosophie DAG

Un DAG (graphe orienté acyclique) permet de structurer un workflow scientifique de manière rigoureuse.

Il apporte plusieurs avantages :

  • Reproductibilité : chaque étape est explicitement définie
  • Performance : mise en cache des calculs intermédiaires
  • Modularité : séparation claire des responsabilités
  • Scalabilité : possibilité d’exécution parallèle

7.9 🔗 Nœud 1 : imports — Initialisation de l’environnement

Rôle dans le DAG : Nœud racine de configuration logicielle

7.9.1 📌 Description

Cette première sous-étape initialise l’ensemble de l’environnement de calcul scientifique nécessaire à l’exécution du pipeline FEM. Elle assure la disponibilité des bibliothèques de calcul numérique, de simulation et de post-traitement.

7.9.2 ⚙️ Fonction

  • Assurer la reproductibilité des simulations
  • Séparer données et implémentation numérique
Code 
# Importation des bibliothèques nécessaires
import numpy as np
import matplotlib.pyplot as plt
from dolfinx import fem, mesh, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import json
from petsc4py import PETSc

print("Bibliothèques importées avec succès!")

import json

def load_problem():
    with open('problem.json', 'r') as f:
        return json.load(f)

problem = load_problem()

# Helper pour convertir les valeurs du JSON (gère "" et null)
def to_float(val):
    if val is None or val == "" or val == "null":
        return None
    return float(val)

def to_int(val):
    if val is None or val == "" or val == "null":
        return None
    return int(val)
Bibliothèques importées avec succès!

7.10 🔗 Nœud 2 : Inputs Parameter — Données du problème

Rôle dans le DAG : Nœud d’ingestion et structuration des données d’entrée

7.10.1 📌 Description

Cette étape est responsable de la lecture, validation et structuration des paramètres du problème physique. Elle transforme une entrée brute (fichier JSON) en structure exploitable par le solveur FEM.

7.10.2 ⚙️ Fonction

  • Définir le cadre physique de simulation
  • Garantir la cohérence des données d’entré
Code 
# --- Données générales ---
L = to_float(problem['length'])
if L is None:
    raise ValueError("L (length) is missing in JSON")

E = problem['material']['E']
if E is None:
    raise ValueError("E (Young modulus) is missing in JSON")
E = float(E)
nu = to_float(problem['material']['nu'])

# --- Supports ---
supportA = problem['supports']['A']
supportB = problem['supports']['B']

# --- Section ---
section = problem["section"]

I = to_float(section["I"])
if I is None:
    raise ValueError("I (moment of inertia) is missing in JSON")

W = to_float(section.get("W"))
h_equiv = to_float(section.get("h_equiv"))
b = to_float(section.get("b"))
h = to_float(section.get("h"))
r = to_float(section.get("r"))
r_int = to_float(section.get("r_int"))
tw = to_float(section.get("tw"))
tf = to_float(section.get("tf"))

# --- Load ---
load = problem["load"]
load_type = load["type"]

P = to_float(load.get("P"))
a = to_float(load.get("xP"))
q = to_float(load.get("q"))
x_start = to_float(load.get("x_start"))
x_end = to_float(load.get("x_end"))

💡 Point clé :

Ces paramètres peuvent être pilotés automatiquement par Link En particulier :

  • 𝑛𝑥,𝑛𝑦,𝑛𝑧 → variables d’optimisation
  • propriétés physiques → constantes

7.11 🔗 Nœud 3 : Mesh — Génération du maillage

Rôle dans le DAG : Nœud d’hyperparamétrisation du modèle numérique

7.11.1 📌 Description

Les paramètres de maillage (nx, ny, nz) ne sont plus considérés comme fixes, mais comme des hyperparamètres du pipeline de simulation.

7.11.2 ⚙️ Fonction

Cette approche permet à l’environnement Link de :

  • explorer automatiquement différentes configurations de maillage
  • évaluer l’impact de la discrétisation sur la précision numérique
  • comparer les solutions obtenues pour différents niveaux de raffinement
  • identifier un compromis optimal entre coût de calcul et précision
Code 
# =========================
# MAILLAGE 3D CENTRÉ
# =========================
#nx = 21
#ny = 11
#nz = ny

msh = mesh.create_box(
    MPI.COMM_WORLD,
    points=[
        np.array([
            0.0,
            -h_equiv / 2,
            -h_equiv / 2
        ]),
        np.array([
            L,
            h_equiv / 2,
            h_equiv / 2
        ])
    ],
    n=[nx, ny, nz],
    cell_type=mesh.CellType.hexahedron
)

7.12 🔗 Nœud 4 : FunctionSpace — Définition de l’espace fonctionnel

Rôle dans le DAG : Nœud de discrétisation des inconnues

7.12.1 📌 Description

Cette étape consiste à définir l’espace fonctionnel dans lequel la solution du problème sera recherchée. Dans le cadre de l’élasticité linéaire 3D, le champ inconnu est le déplacement vectoriel. On introduit donc un espace de fonctions vectorielles basé sur des éléments finis de type Lagrange. Les fonctions de test et d’essai sont également définies à ce stade, conformément à la formulation variationnelle.

7.12.2 ⚙️ Fonction

  • Définir les degrés de liberté du problème
  • Introduire les inconnues et fonctions de test
  • Préparer la formulation faible

7.13 🔗 Nœud 5 : CLs — Définition de l’espace fonctionnel

Rôle dans le DAG : Nœud de contrainte du système

7.13.1 📌 Description

Les conditions aux limites sont imposées afin de garantir l’unicité et la pertinence physique de la solution. Dans ce cas, un encastrement est appliqué à une extrémité de la poutre, imposant un déplacement nul.

7.13.2 ⚙️ Fonction

  • Contraindre le système mécanique
  • Représenter les appuis physiques
  • Assurer la stabilité numérique
Code 
# =========================
# CONDITIONS AUX LIMITES (encastrement x=0)
# =========================
def left(x):
    return np.isclose(x[0], 0.0)

left_dofs = fem.locate_dofs_geometrical(V, left)

bc = fem.dirichletbc(
    np.array((0.0, 0.0, 0.0), dtype=PETSc.ScalarType),
    left_dofs,
    V
)

7.14 🔗 Nœud 6 : Loads — Application des chargements

Rôle dans le DAG : Nœud de sollicitation externe

7.14.1 📌 Description

Cette étape introduit les actions mécaniques appliquées à la structure. Dans ce cas, une charge est appliquée localement sur une zone de la surface. La force est convertie en traction surfacique équivalente, compatible avec la formulation faible.

7.14.2 ⚙️ Fonction

  • Modéliser les sollicitations physiques
  • Injecter l’énergie dans le système
  • Définir le second membre du problème
Code 
# =========================
# CHARGEMENTS PONCTUELLE
# =========================

tol = L / ny

def load_patch(x):
    return (
        (np.abs(x[0] - L) < tol) &
        (np.abs(x[1]) <= 0.51 * h_equiv) &
        (x[2] > 0.49 * h_equiv)
    )

fdim = msh.topology.dim - 1

facets = mesh.locate_entities_boundary(msh, fdim, load_patch)

facet_tag = mesh.meshtags(
    msh,
    fdim,
    facets,
    np.full(len(facets), 1, dtype=np.int32)
)

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)

area = fem.assemble_scalar(fem.form(1 * ds(1)))

g = fem.Constant(
    msh,
    PETSc.ScalarType((0.0, 0.0, -P / area))
)

7.15 🔗 Nœud 7 : MechanicalBehaviour — Comportement mécanique

Rôle dans le DAG : Nœud de modélisation constitutive du matériau

7.15.1 📌 Description

Ce nœud définit la loi de comportement du matériau, c’est-à-dire la relation entre les déformations et les contraintes dans la structure. La relation contrainte–déformation est donnée par la loi de Hooke généralisée, formulée à l’aide des coefficients de Lamé.

7.15.2 ⚙️ Fonction

  • Définir la physique du matériau
  • Relier le champ de déplacement aux contraintes internes
  • Fournir l’expression du tenseur des contraintes pour la formulation faible
Code 
# =========================
# MÉCANIQUE
# =========================
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u):
    return (
        lmbda * ufl.tr(epsilon(u)) * ufl.Identity(3)
        + 2 * mu * epsilon(u)
    )

7.16 🔗 Nœud 8 : Formumation — Formulation faible

Rôle dans le DAG : Nœud de traduction mathématique du problème physique

7.16.1 📌 Description

Ce nœud définit la loi de comportement du matériau, c’est-à-dire la relation entre les déformations et les contraintes dans la structure. La relation contrainte–déformation est donnée par la loi de Hooke généralisée, formulée à l’aide des coefficients de Lamé.

7.16.2 ⚙️ Fonction

  • Définir la physique du matériau
  • Relier le champ de déplacement aux contraintes internes
  • Fournir l’expression du tenseur des contraintes pour la formulation faible
Code 
# =========================
# FORMULATION FAIBLE
# =========================
a_form = ufl.inner(sigma(u), epsilon(v)) * ufl.dx(domain=msh)

f = fem.Constant(
    msh,
    PETSc.ScalarType((0.0, 0.0, 0.0))
)

L_form = ufl.inner(g, v) * ds(1)

7.17 🔗 Nœud 9 : Solve — Résolution du système FEM

Rôle dans le DAG : Nœud principal de calcul

7.17.1 📌 Description

Le problème variationnel est discrétisé et transformé en un système linéaire de la forme \(𝐾𝑈=𝐹\). Ce système est ensuite résolu à l’aide d’un solveur basé sur PETSc.

7.17.2 ⚙️ Fonction

  • Assembler le système matriciel
  • Appliquer les conditions aux limites
  • Calculer le champ de déplacement
Code 
# =========================
# RÉSOLUTION
# =========================
problema = LinearProblem(
    a_form,
    L_form,
    bcs=[bc],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu"
    }
)

uh = problema.solve()

print("✅ Problème résolu")

# =========================
# POST SAFE
# =========================
vals = uh.x.array
w_max = np.min(vals)
print("Flèche max (mm) =", w_max * 1000)
✅ Problème résolu
Flèche max (mm) = -0.011385892761204284

7.18 🔗 Nœud 10 : PostProcess — Post-traitement mécanique

Rôle dans le DAG : Nœud d’analyse des champs mécaniques

7.18.1 📌 Description

À partir du champ de déplacement obtenu, on calcule des grandeurs physiques dérivées, en particulier la contrainte équivalente de von Mises.

7.18.2 ⚙️ Fonction

  • Évaluer l’état de contrainte dans la structure
  • Identifier les zones critiques
  • Préparer les résultats pour visualisation
Code 
# =========================
# POST-TRAITEMENT
# =========================
# Von Mises
sigma_u = sigma(uh)

dev = sigma_u - (1/3) * ufl.tr(sigma_u) * ufl.Identity(3)
von_Mises = ufl.sqrt(3/2 * ufl.inner(dev, dev))

Vsig = fem.functionspace(msh, ("Lagrange", 1))
vm = fem.Function(Vsig)

vm_expr = fem.Expression(von_Mises, Vsig.element.interpolation_points())
vm.interpolate(vm_expr)

# =========================
# EXPORT PARAVIEW
# =========================
#with io.XDMFFile(MPI.COMM_WORLD, "beam_results.xdmf", "w") as file:
#    file.write_mesh(msh)
#    uh.name = "Displacement"
#    vm.name = "Von_Mises"
#    file.write_function(uh)
#    file.write_function(vm)

#print("✅ Résultats exportés (Paraview)")

7.19 🔗 Nœud 11 : Visualisation — Visualisation des réultats

Rôle dans le DAG : Nœud de restitution et d’interprétation visuelle

7.19.1 📌 Description

Les résultats sont projetés sur une représentation graphique du maillage afin de visualiser :

  • la déformée de la structure
  • le champ de déplacement
  • le champ de contraintes

7.19.2 ⚙️ Fonction

  • Faciliter l’interprétation des résultats
  • Valider qualitativement le comportement mécanique
  • Produire des supports visuels pour l’analyse

7.20 🔗 Nœud 12 : Comaparison — Comparaison à la théorie

Rôle dans le DAG : Nœud d’analyse ciblée et de validation physique

7.20.1 📌 Description

Ce nœud calcule l’écart entre la solution numérique obtenue par la méthode des éléments finis et une solution de référence analytique. Il introduit une métrique quantitative de précision, appelée erreur relative, qui permet d’évaluer la qualité de la simulation.

7.20.2 ⚙️ Fonction

  • Valider localement la solution numérique
  • Comparer FEM 3D vs modèle simplifié (RDM)
  • Mettre en évidence les écarts liés au maillage

FEniCS-Link

7.21 🔗 Nœud 13 : error — Évaluation de l’erreur

Rôle dans le DAG : Nœud d’évaluation et de pilotage de l’optimisation

7.21.1 📌 Description

Ce nœud calcule l’écart entre la solution numérique obtenue par la méthode des éléments finis et une solution de référence analytique. Il introduit une métrique quantitative de précision, appelée erreur relative, qui permet d’évaluer la qualité de la simulation.

7.21.2 ⚙️ Fonction

  • Quantifier la précision du modèle FEM
  • Comparer différentes configurations de maillage
  • Fournir un critère objectif pour l’optimisation

L’erreur relative est définie par :

\[ \varepsilon = \frac{|y_{\max}^{FEM} - y_{\max}^{ref}|}{y_{\max}^{ref}} \]

Code 
# =========================
# COMPARAISON
# =========================

y_fem_at_L = w_max*1000
y_analytique_at_L=-0.01193459541113796

erreur_relative = abs(y_fem_at_L - y_analytique_at_L) / abs(y_analytique_at_L) * 100

# Resultats de comparaison
print("=== Comparaison flèche max. ===")
print(f"Analytique: {y_analytique_at_L*1000:.5f} mm")
print(f"FEM 3D:     {y_fem_at_L*1000:.5f} mm")
print(f"Erreur:     {erreur_relative:.2f} %")

# target value (HPO):
erreur_relative
=== Comparaison flèche max. ===
Analytique: -11.93460 mm
FEM 3D:     -11.38589 mm
Erreur:     4.60 %
np.float64(4.597580655491677)

7.22 🧠 Résumé du projet

Ce travail s’inscrit dans l’exploration d’un environnement de calcul scientifique basé sur des graphes orientés acycliques (DAG) appliqués à la simulation par éléments finis et à l’optimisation automatique de paramètres.

L’objectif principal est d’utiliser une approche de type workflow pour structurer un pipeline complet allant de la définition du problème mécanique jusqu’à l’optimisation des paramètres de maillage.

7.22.1 ⚙️ Environnement de calcul basé sur DAG

L’extension Link pour JupyterLab permet de transformer un notebook classique en un graphe de calcul structuré, dans lequel chaque cellule devient un nœud du DAG.

Cette approche apporte plusieurs avantages :

  • exécution contrôlée des dépendances
  • mise en cache des résultats intermédiaires
  • réduction des recomputations inutiles
  • meilleure traçabilité des pipelines numériques
  • possibilité de parallélisation de certaines étapes

Dans ce cadre, un workflow scientifique n’est plus linéaire, mais devient un système structuré et reproductible.

7.22.2 🔬 Application au problème FEM

Le cas étudié est celui d’une poutre en flexion résolue par éléments finis (FEniCSx).

Le pipeline complet inclut :

  • génération du maillage paramétré par \((n_x, n_y, n_z)\)
  • résolution du problème mécanique linéaire
  • calcul du champ de déplacement
  • post-traitement (flèche, contraintes de von Mises)
  • extraction de métriques d’erreur

7.22.3 📊 Optimisation et étude de convergence

Les paramètres de maillage \((n_x, n_y, n_z)\) sont considérés comme des hyperparamètres.

La sortie du système est une fonction objectif :

  • erreur relative entre solution FEM et solution analytique

Cette erreur est utilisée dans un processus de Hyperparameter Optimization (HPO) afin de :

  • minimiser l’erreur numérique
  • analyser la convergence du maillage
  • identifier un compromis précision / coût de calcul

7.22.4 🧩 Apport du DAG dans ce contexte

L’utilisation d’un DAG permet de structurer naturellement le workflow :

  • chaque étape de la simulation devient un nœud
  • les dépendances sont explicites
  • les résultats sont réutilisables via caching
  • les modifications locales ne déclenchent pas tout le pipeline

Cela transforme une étude FEM classique en un pipeline computationnel reproductible et optimisable.

7.23 💡 Conclusion

Ce projet illustre comment un environnement basé sur DAG peut améliorer les workflows de simulation numérique en :

  • automatisant les études de convergence
  • structurant les calculs FEM
  • facilitant l’intégration de méthodes d’optimisation (HPO)
  • réduisant les coûts de recalcul

Il ouvre ainsi la voie vers des approches plus robustes de simulation scientifique pilotée par graphes et optimisation automatique.

7.23.1 🐙 Accès au projet et apprentissage

L’ensemble du projet est disponible en open source sur GitHub : fenics-link et peut être exécuté directement via un environnement Docker, sans installation locale.

Les utilisateurs souhaitant tester concrètement le workflow (DAG + FEM + optimisation) peuvent ainsi explorer l’outil de manière interactive et intuitive.

Pour des raisons de lisibilité du document, les démonstrations détaillées de l’interface Link et des expérimentations HPO ne sont pas reproduites ici.

7.23.2 🎓 Ressources et accompagnement

Pour les utilisateurs souhaitant une démonstration guidée ou une prise en main plus approfondie :

  • une page de contact est disponible pour obtenir des explications complémentaires : contact
  • des tutoriels et webinaires peuvent être proposés sur demande
  • un webinaire sera organisé dès qu’un seuil d’intérêt suffisant sera atteint (≈ 10 participants)