.. _sec:opti_topo_intro:

Éléments théoriques de base
===========================

Généralités
-----------

L'optimisation topologique consiste à trouver la **répartition optimale de la matière** dans une **pièce donnée**
soumise à des **chargements** et sous certaines **limitations**. L'inconnue de ce type de problème est la *topologie*,
c'est-à-dire *comment la matière est-elle répartie ?*

Dans un cadre mécanique, ce type de problème peut être formulé en : **trouver la distribution optimale de rigidité**
:math:`\mathbfcal{k}`, ce qui peut s'écrire :

.. math::
   :name: eq:opti_topo_1

   \min_{\mathbfcal{k}(m)} \quad & f(\mathbfcal{k})   & \\
   \textsf{tel que}        \quad & g_i(\mathbfcal{k}) &= 0 \\

Explicitons certains termes :

- Par **répartition optimale de la matière** on entend la distribution spatiale de rigidité :math:`\mathbfcal{k}(m)`
  en tout point :math:`m \in \Omega` qui minimise la *fonction objectif* :math:`f`, indicatrice du souhait de l'utilisateur.
- Par **pièce donnée** on entend que la recherche de la topologie se fait sur un domaine limité de l'espace :math:`\Omega`.
- Par **chargements** on entend que cette pièce subit des conditions aux limites, par exemple des déplacements
  imposés et/ou des forces appliquées sur certaines zones.
- Par **limitations** on entend que l'optimisation se fait *sous contraintes*, représentées par l'ensemble de fonctions
  :math:`g_i(\mathbfcal{k}) = 0`.

Le problème est illustré ci-dessous dans le cas d'une poutre en flexion 3 points.

.. figure:: figures/99_lines_img_1.png
   :name: fig:opti_topo_flexion
   :width: 15cm
   :align: center

   Problème d'optimisation topologique sur une poutre en flexion [SIGMUND-2001]_

Plusieurs choix de fonctions objectif :math:`f` et contraintes :math:`g_i` sont alors possibles pour optimiser
une structure, comme par exemple :

- chercher la pièce la plus rigide possible, sans dépasser un certain volume ;
- chercher la pièce la plus légère possible, sans dépasser une certaine fréquence de résonnance ;
- chercher la pièce avec les contraintes mécaniques les plus faibles, sans dépasser un certain niveau de déplacement.

Un choix très répendu est de **minimiser la compliance** (c'est-à-dire *l'énergie de déformation* ou encore *le travail
des forces extérieures*) **sous contrainte de volume**.

Dans ce document nous illustrerons brièvement la méthode **SIMP** (pour *Solid Isotropic Material with Penalization*)
qui est la méthode d'optimisation topologique la plus répendue dans les codes industriels et celle mise en oeuvre
dans Cast3M via la procédure `TOPOPTIM <http://www-cast3m.cea.fr/index.php?page=notices&notice=TOPOPTIM>`_.
Le lecteur intéressé pourra consulter de nombreux ouvrages sur le sujet, comme par exemple :

- Les livres de référence [BENDSOE-1995]_ et [BENDSOE-AL-2004]_ qui détaillent rigouresement la théorie
  derrière l'optimisation topologique.
- Les articles pédagogiques [SIGMUND-2001]_ et [ANDREASSEN-AL-2011]_ qui présentent des implémentations sur Matlab,
  en respectivement 99 lignes puis 88 lignes, d'un algorithme d'optimisation topologique. La procédure
  `TOPOPTIM <http://www-cast3m.cea.fr/index.php?page=notices&notice=TOPOPTIM>`_ de Cast3M en est inspiré et
  l'exemple utilisé dans ce document en est une reproduction.

.. _sec:opti_topo_simp:

La méthode SIMP
---------------

Le domaine de conception est discrétisé selon un maillage de :math:`N` éléments. Chaque élément :math:`e` se voit attribué
une *variable de densité* :math:`x_e` et un module d'Young :math:`E_e` selon une loi puissance :

.. math::
   :name: eq:loi_simp

   E_e(x_e) = E_{\textrm{min}} + {x_e}^p(E_0-E_{\textrm{min}})

avec :

- :math:`x_e \in [0,1]` la densité de l'élément :math:`e` (où :math:`x_e=0` correspond à un élément "vide" et :math:`x_e=1` un élément "plein")
- :math:`E_0` le module d'Young du matériau
- :math:`E_{\textrm{min}}` une valeur "faible" de module d'Young choisie pour éviter que la matrice de rigidité du modèle éléments finis devienne singulière (typiquement :math:`E_{\textrm{min}}=E_0.10^{-6}`)
- :math:`p` un coefficient de *pénalisation* (typiquement :math:`p=3`) introduit afin de défavoriser les valeurs intermédiaires :math:`0<x_e<1` lors de l'optimisation

Le problème d'optimisation topologique de **minimisation de la compliance**, sous contrainte de volume, s'écrit alors :

.. math::
   :name: eq:opti_topo_2

   \min_{\textbf{x}} \quad & f(\textbf{x}) = \textbf{U}^T.\mathbfcal{K}.\textbf{U}
                                           = \sum_{e=1}^{N} E_e(x_e)\textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e \\
   \textsf{tel que}  \quad & g(\textbf{x}) = \phi(\textbf{x}) - \phi_0 = \sum_{e=1}^N \frac{x_eV_e}{V_0} - \phi_0 = 0 \\
                           & \mathbfcal{K}.\textbf{U} = \textbf{F} \\
                           & 0 \leq x_e \leq 1 \quad \forall e \in \{1;\dots;N\} \\

Ce problème de minimisation de la compliance est équivalent à *maximiser la raideur de la structure*. Précisons quelques variables :

- :math:`\textbf{x}` est le vecteur des densités des :math:`N` éléments du maillage, les iconnues du problème, encore appellées *variables de conception*
- :math:`f` est la fonction objectif, la *compliance* de la structure (homogène à une énergie)
- :math:`\textbf{U}` et :math:`\textbf{F}` sont les vecteurs déplacements et forces aux noeuds du maillage
- :math:`\textbf{u}_e` est le vecteur déplacements réduit à l'élément :math:`e`
- :math:`\mathbfcal{K}=\sum_{e=1}^N E(x_e)\mathbfcal{k}_0` est la matrice de rigidité globale, assemblée sur les :math:`N` éléments du maillage
- :math:`\mathbfcal{k}_0` est la matrice de rigidité élémentaire pour un module d'Young unitaire
- :math:`g(\textbf{x})` est la fonction contrainte d'égalité
- :math:`\phi(\textbf{x})` est la fraction volumique de la topologie et :math:`\phi_0` est la fraction volumique imposée
- :math:`V_e` est le volume de l'élément :math:`e` et :math:`V_0` est le volume du domaine de conception :math:`\Omega`

.. _sec:opti_topo_oc:

Résolution par Critère d'Optimalité
-----------------------------------

Un schéma de résolution heuristique et simple de ce type du problème :eq:`eq:opti_topo_2` est
proposé par [BENDSOE-1995]_ et consite à mettre à jour, de manière itérative, les densités
courantes :math:`\textbf{x}` vers une nouvelle valeur :math:`\textbf{x}^{\textrm{new}}` :

.. math::
   :name: eq:opti_topo_bendsoe

   x_e^{\textrm{new}} = \left\{
     \begin{array}{lll}
     x_e^-         & \textsf{si} & x_eB_e^{\eta} \le x_e^- \\
     x_eB_e^{\eta} & \textsf{si} & x_e^- < x_eB_e^{\eta} < x_e^+ \\
     x_e^+         & \textsf{si} & x_e^+ \le x_eB_e^{\eta} \\
     \end{array}
   \right.

où :

- :math:`x_e^- = \max (0,x_e-m)` est une borne inférieure pour la densité à l'itération courante
- :math:`x_e^+ = \min (1,x_e+m)` est une borne supérieure pour la densité à l'itération courante
- :math:`m` est la limite d'incrément de densité à chaque itération pour stabiliser la convergence
- :math:`\eta` est un coefficient d'amortissement (généralement :math:`\eta=0,5`)

Le terme :math:`B_e` guidant la mise à jour de :math:`x_e` est obtenu par la condition d'optimalité :

.. math::
   :name: eq:opti_topo_optimalite

   B_e = \frac{-\dfrac{\partial f}{\partial x_e}}{\mathcal{L} \dfrac{\partial g}{\partial x_e}}

- :math:`\dfrac{\partial f}{\partial x_e}` est la **sensibilité** de la fonction objectif :math:`f`
- :math:`\dfrac{\partial g}{\partial x_e}` est la **sensibilité** de la fonction contrainte :math:`g`
- :math:`\mathcal{L}` est le **multiplicateur de Lagrange** à déterminer pour satisfaire la contrainte :math:`g`

En dérivant les expressions des fonctions, la sensibilité de la fonction objectif (compliance), en supposant l'absence
de forces dépendantes de la densité :math:`(\frac{\partial \textbf{F}}{\partial x_e}=0)`, s'écrit :

.. math::
   :name: eq:opti_topo_sensibilite_1

   \frac{\partial f}{\partial x_e} = -p(x_e)^{p-1}(E_0-E_{\textrm{min}}) \textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e

La sensibilité de la fonction contrainte (fraction volumique) s'écrit :

.. math::
   :name: eq:opti_topo_sensibilite_2

   \frac{\partial g}{\partial x_e} = \frac{V_e}{V_0}

La difficulté étant alors de trouver la valeur de :math:`\mathcal{L}` qui satisfait la contrainte :math:`g`. Celle-ci
peut être calculée par un algorithme de **dichotomie** en initialisant des bornes inférieure :math:`\mathcal{L}^-` et
supérieure :math:`\mathcal{L}^+` puis en choisissant la valeur milieu de l'intervalle. Une évaluation de la fonction
contrainte :math:`g` est alors faite et le processus est répété dans le demi intervalle *ad hoc* :

.. _algo:opti_topo_dichotomie:

**Initialisation des bornes**

:math:`\mathcal{L}^- =0 \quad \mathcal{L}^+ =10^9`

**Tant que** \ :math:`\frac{\mathcal{L}^+ - \mathcal{L}^-}{\mathcal{L}^+ + \mathcal{L}^-} > 10^{-3}` :

.. raw:: html

   <div style="margin-left:20px;width:300px;height:240px;">
.. math::

   \begin{array}{ll}
     \mathcal{L}                   & = (\mathcal{L}^- + \mathcal{L}^+)/2 \\
     \textbf{B}                    & = \textsf{calculer selon (5)} \\
     \textbf{x}^{\textrm{new}}     & = \textsf{calculer selon (4)} \\
     g(\textbf{x}^{\textrm{new}})  & = \textsf{calculer selon (3)} \\
     \textsf{si } g(\textbf{x}^{\textrm{new}}) & > 0 & \\
       \quad \mathcal{L}^- & = \mathcal{L} \\
     \textsf{sinon} & \\
       \quad \mathcal{L}^+ & = \mathcal{L} \\
     \textsf{finsi} &\\
   \end{array} \\

.. raw:: html

   </div>

**Fin**

À l'issue de la dichotomie on obtient la valeur de :math:`\mathcal{L}` qui satisfait la contrainte
sur la fraction volumique ainsi que la nouvelle topologie :math:`\textbf{x}^{\textrm{new}}`.

.. _sec:opti_topo_filtre_sens:

Filtrage de la sensibilité
--------------------------

Afin d'éviter l'effet de damier et diminuer la sensibilité des solutions au maillage, on applique une
procédure de filtrage (ou lissage) sur le champ de desnité ou sur les champs de sensibilité.

Le **filtrage de la sensibilité** consiste à modifier la sensibilité :math:`\partial f / \partial x_e` d'un
élément :math:`e` par une valeur moyenne pondérée des sensibilités des éléments voisins :

.. math::
   :name: eq:opti_topo_filtre_sens

   \dfrac{\widehat{\partial f}}{\partial x_e} = \frac{1}{\textrm{max}(10^{-3},x_e)}\frac{\sum_{i \in N_e}\hat{H}_{ei}x_i\frac{\partial f}{\partial x_i}}{\sum_{i \in N_e}\hat{H}_{ei}}

Les coefficients :math:`\hat{H}_{ei}` valent :

.. math::

   \hat{H}_{ei} = \left( 1 - \frac{\textrm{dist}(e,i)}{r_{\textrm{min}}} \right)^q V_i

et n'est définit que pour les :math:`N_e` éléments :math:`i` tels que :math:`\textrm{dist}(e,i) \le r_{\textrm{min}}`.
Détaillons les nouvelles variables :

- :math:`\textrm{dist}(e,i)` la distance entre les centres des éléments :math:`e` et :math:`i`
- :math:`V_i` le volume de l'élément :math:`i` (ou bien une autre quantitié pour pondérer)
- :math:`r_{\textrm{min}}` le rayon du filtre, au dela duquel l'opérateur :math:`\hat{H}_{ei}` est nul
- :math:`q` un coefficient

Dans Cast3M, ce filtrage est réalisé grâce aux opérateurs `MFIL <http://www-cast3m.cea.fr/index.php?page=notices&notice=MFIL>`_
(pour calculer le terme au numérateur) et `NFIL <http://www-cast3m.cea.fr/index.php?page=notices&notice=NFIL>`_ (pour normaliser
par le terme au dénominateur).

Notons que dans les articles des codes 99 lignes [SIGMUND-2001]_ et 88 lignes [ANDREASSEN-AL-2011]_, l'opérateur de filtrage
utilisé correspond au cas où :math:`q=1` et où tous les éléments ont un volume :math:`V_i=1`.

.. _sec:opti_topo_filtre_dens:

Filtrage des densités
---------------------

*A faire ...*

.. _sec:opti_topo_seuil:

Seuillage des densités
----------------------

*A faire ...*