\documentclass[11pt]{amsart}
\usepackage{amsfonts}
%\usepackage{french}
\usepackage{amscd}
\usepackage[dvips]{epsfig}

\pagestyle{empty}
\oddsidemargin 0.8cm
\evensidemargin 0.2cm
\textwidth 15.2 cm
\topmargin 0cm
\textheight 21cm

\begin{document}
\author{Philippe Guillaume, Mohamed Masmoudi}
\title{TP Conception Optimale de Forme}
\date{}
\maketitle
 
\centerline{DEA de Math\'ematiques Appliqu\'ees} 
\centerline{Ecole Doctorale de Math\'ematiques de Toulouse} 
\centerline{INSA - SUPAERO - UPS - UT1} 

\vspace{2cm} 
\begin{figure}[h]
\centerline{\ \psfig{file=poutre3.eps,width=10cm,height=6cm} } %\caption{}
\end{figure}

\newpage

\section{D\'{e}finition du TP}

Le but de ce TP est d'optimiser la forme de quelques structures 2D soumises
\`{a} un chargement. Le crit\`{e}re \`{a} optimiser est l'\'energie de
d\'eformation. Pour cela, partant d'un rectangle plein, on utilise le
gradient topologique pour retirer de la mati\`{e}re progressivement jusqu'\`{a}
une proportion fix\'{e}e par l'utilisateur. Nous \'{e}tudierons trois cas :

\begin{enumerate}
\item  le cas d'une poutre de longueur 0,4 et de hauteur 1, avec condition
de Dirichlet nulle sur le bord gauche. Une force dirig\'{e}e vers le bas est
appliqu\'{e}e au milieu du c\^{o}t\'{e} oppos\'{e}. La solution de ce 
probl\`{e}me est connue : c'est une \'{e}querre,

\item  le cas d'une poutre de longueur 1,6 et de hauteur 1, avec les
m\^{e}mes conditions aux limites. Les solutions de ce probl\`{e}me sont 
connues : ce sont les treillis de Michell,

\item  le cas d'un pont de longueur 2 et de hauteur 1, avec condition de
Dirichlet sur le coin inf\'{e}rieur droit, et condition de glissement
horizontal sur le coin inf\'{e}rieur gauche. Une force dirig\'{e}e vers le
bas est appliqu\'{e}e au milieu du cot\'{e} inf\'{e}rieur. Une solution
connue est en forme de demie roue avec rayons.
\end{enumerate}

\section{D\'{e}roulement du TP}

Ce TP est effectu\'{e} sous {\sc matlab} avec utilisation de la 
{\sc pde Toolbox}. Les diff\'{e}rents fichiers dont vous aurez besoin sont 
dans le r\'{e}pertoire {\tt /echange/tpcof}. 
Les fichiers {\tt $*$.m} sont des fichiers source {\sc matlab}, les fichiers 
{\tt $*$.mat} sont des fichiers de donn\'{e}es binaires.

Pour chacun des trois cas \'{e}tudi\'{e}s, vous aurez \`{a} :

\begin{itemize}
\item  construire le maillage du rectangle via l'interface de la {\sc pde
Toolbox}. Pour cela, lancez la commande \textsc{matlab}, puis ex\'{e}cutez
l'un des trois fichiers \texttt{rect0, poutre0} ou\texttt{\ pont0.}
L'interface de la {\sc pde toolbox} appara\^{\i }t. La g\'{e}om\'{e}trie, les
propri\'{e}t\'{e}s de mat\'{e}riau, les coefficients de l'EDP et les
conditions aux limites sont pr\'{e}d\'efinies. Choisissez alors un pas de
maillage dans le menu \textsc{mesh/parameters,} et lancez le mailleur avec
le menu \textsc{mesh/initialze mesh} suivi de quelques \textsc{mesh/jiggle
mesh} pour am\'{e}liorer la qualit\'{e} du maillage,

\item  exporter les donn\'{e}es dans l'environnement {\sc matlab}. Pour cela,
choisissez la commande \textsc{export} dans chacun des menus \textsc{draw,
boundary, pde, mesh.} En utilisant la commande \texttt{whos} dans la
fen\^{e}tre {\sc matlab}, vous devez voir appara\^{\i }tre les noms 
(par d\'efaut) des variables 
\textit{a, b, c, d, e, f, g, gd, ns, p, sf, t}. Vous pouvez les sauvegarder
avec la commande \texttt{save fichier. }Pour les r\'{e}cup\'{e}rer
ult\'{e}rieurement, vous utiliserez la commande \texttt{load fichier.}
\end{itemize}

Vous ex\'{e}cuterez ensuite votre programme d'optimisation. Un exemple de
programme (\texttt{main0.m}) se trouve dans le r\'{e}pertoire 
\texttt{/echange/tpcof}. A ne consulter qu'en cas d'extr\^{e}me 
n\'{e}cessit\'{e} !

\section{Description du probl\`{e}me}

\subsection{Formulation continue}

Soit $\Omega $ un domaine born\'{e} de $\Bbb{R}^{2}$. Sa fronti\`{e}re $%
\Gamma $ est compos\'{e}e de deux parties disjointes $\Gamma _{1}$ et $%
\Gamma _{2}$. Le champ $u$ des d\'{e}placements est solution des
\'{e}quations de l'\'{e}lasticit\'{e} lin\'{e}aire 
\begin{equation*}
\left\{ 
\begin{array}{c}
-\mathrm{div\,}\sigma (u)=0\quad \mathrm{dans\ }\Omega , \\ 
u=0\quad \mathrm{sur\ }\Gamma _{1}, \\ 
\sigma (u)n=g\quad \mathrm{sur\ }\Gamma _{2},
\end{array}
\right.
\end{equation*}
avec les notations suivantes : 
\begin{eqnarray*}
\varepsilon (u) &=&\frac{1}{2}(Du+Du^{T}), \\
\sigma (u) &=&hH_{0}\varepsilon (u),
\end{eqnarray*}
o\`{u} $H_{0}$ est le tenseur de Hooke du mat\'{e}riau, constant sur le
domaine, et $h(x)=$ 1 ou $\epsilon ,$ la variable de programmation $\epsilon 
$ \'{e}tant de l'ordre de $10^{-3}.$

Le probl\`{e}me d'optimisation de forme est le suivant :

\begin{eqnarray*}
&&\left\{ 
\begin{array}{l}
\min \ J(h) \\ 
\int_{\Omega }h(x)\,dx=C,
\end{array}
\right. \\
J(h) &=&\int_{\Omega }\sigma (u):\varepsilon (u)\,dx=\int_{\Gamma
_{2}}g.u\,dx,
\end{eqnarray*}
o\`{u} $C$ est la proportion de mat\'{e}riau \`{a} conserver, fix\'{e}e par
vous.

Le gradient topologique est approch\'e en chaque point par 
\begin{equation*}
J_{T}^{\prime }(x)=2\,\sigma (u):\varepsilon (u).
\end{equation*}
Une it\'{e}ration de l'algorithme d'optimisation consiste \`{a} retirer une
partie $P$ du domaine telle que 
\begin{equation*}
J_{T}^{\prime }(x)\leq J_{T}^{\prime }(y),\quad \forall (x,y)\in P\times
\Omega \backslash P.
\end{equation*}
La proportion mes($P$)/mes($\Omega $) est une variable de programmation
\`{a} fixer entre 0 et 0.3. On r\'{e}initialise alors $\Omega $ par $\Omega
\leftarrow \Omega \backslash P$ et $h$ par $h(x)=\epsilon $ si $x\in P,$ $%
h(x)=1$ sinon. L'algorithme s'arr\^{e}te lorsque la proportion voulue de
mat\'{e}riau est atteinte, la solution obtenue \'{e}tant le domaine $\Omega
\backslash P.$

\subsection{Formulation discr\`{e}te}

On utilise une m\'{e}thode $P_{1}$ sur un maillage non structur\'{e}. Le
d\'{e}placement discret $u$ est continu et affine par morceaux, et on note $%
U $ le vecteur des coordonn\'{e}es de $u$ dans la base des \'{e}l\'{e}ments
finis.

La matrice de rigidit\'{e} est not\'{e}e $K,$ elle d\'{e}pend de $H_{0}$
(constant sur tout le domaine) et de $h$ (constant sur chaque
\'{e}l\'{e}ment). Le chargement $g$ d\'{e}finit le second membre $G$, et le
vecteur $U=U(h)$ est solution du syst\`{e}me 
\begin{equation*}
KU=G.
\end{equation*}
Le probl\`{e}me d'optimisation discret devient 
\begin{equation*}
\left\{ 
\begin{array}{l}
\min \ G.U \\ 
\int_{\Omega }h(x)\,dx=C.
\end{array}
\right.
\end{equation*}

\subsection{Quelques commandes utiles de {\sc matlab}}

\begin{itemize}
\item  {\tt help, help help, help topic,... }: aide en ligne,

\item  {\tt help sparfun } : aide sur les matrices creuses,

\item  {\tt whos } : affiche vos variables courantes,

\item  {\tt a.$*$b } : multiplication terme \`{a} terme des vecteurs ou des
matrices $a$ et $b$,

\item  {\tt M' } :\ transpos\'{e}e (conjugu\'{e}e) de la matrice M,

\item  {\tt [nl,nc] = size(M) } : nombre de lignes et de colonnes de la matrice M,

\item  {\tt programme } : si vous avez cr\'{e}\'{e} un fichier nomm\'{e}
{\tt programme.m} et situ\'{e} dans votre r\'{e}pertoire courant, la commande
{\tt programme} ex\'{e}cute les lignes de votre programme,

\item  {\tt [x,y,...] = fonction(a,b,...) } : si vous avez cr\'{e}\'{e} 
dans votre r\'{e}pertoire courant un fichier 
nomm\'e {\tt fonction.m} ayant pour ent\^{e}te\\
\null \qquad {\tt function [x,y,...] = fonction(a,b,...)}\\
la commande {\tt [x,y,...] = fonction(a,b,...)} vous renvoie les valeurs
{\tt [x,y,...]} calcul\'{e}es par votre fonction \`{a} partir des valeurs
d'entr\'{e}e {\tt (a,b,...)}. Les variables {\tt x,y,a,b,...} peuvent 
\^{e}tre de tailles diff\'{e}rentes,

\item  {\tt [y,I] = sort(x)} : ordonne par ordre croissant les \'{e}l\'{e}ments
du vecteur {\tt x} dans le vecteur {\tt y}. Les indices {\tt I} sont tels que 
{\tt y = x(I)}.

\item  {\tt disp(x)} : affiche la valeur de la variable $x$,

\item  {\tt disp('cha\^{\i }ne')} : affiche la cha\^{\i }ne de caract\`{e}res
{\tt cha\^{\i }ne},

\item  {\tt pdeplot(p,e,t,'xydata',ep,'xystyle','flat','mesh','on')} :
visualisation sur le maillage des donn\'{e}es contenues dans le vecteur {\tt ep}
(cf. {\tt help pdeplot}),

\item  {\tt pause(1)} : arr\^{e}t buffet,

\item  les boucles :\newline
{\tt for i = 1:10

instructions;\newline
end\newline
while val $>=$ 0

instructions;\newline
end}
\end{itemize}

\subsection{Concernant la r\'{e}solution de l'edp}

\begin{itemize}
\item  Sous {\sc matlab}, le tenseur $\sigma (u)$ s'\'{e}crit $c\otimes \nabla u.$
La variable $c$ que vous avez export\'{e}e au d\'{e}part est une cha\^{\i %
}ne de caract\`{e}res, que vous transformerez en une matrice de dimension
{\tt 10 $\times $ nbt} ({\tt nbt} = nombre de triangles) par :\newline
{\tt c0 = zeros(1,10);\newline
for i = 1:10

c0(i) = eval(c(i,:));\newline
end\newline
c = diag(c0)$*$ones(10,nbt);}\newline
Chaque colonne de $c$ contient alors les informations relatives au tenseur
de Hooke $H_{0}$ sur l'\'{e}l\'{e}ment correspondant, \`{a} multiplier
ensuite par $h(x).$

\item  {\tt [K,M,F,Q,G,H,R] = assempde(b,p,e,t,c,a,f)} : calcul de la matrice de
rigidit\'{e} $K$ et du second membre $G.$

\item  {\tt U = assempde(K,M,F,Q,G,H,R)} : calcul de la solution.

\item  {\tt [ux,uy] = pdegrad(p,t,U)} : sur le triangle num\'{e}ro $j,$ on a 
\begin{equation*}
Du(x)=\left[ 
\begin{array}{cc}
\mathrm{ux}(1,j) & \mathrm{uy}(1,j) \\ 
\mathrm{ux}(2,j) & \mathrm{uy}(2,j)
\end{array}
\right] ,
\end{equation*}

\item  {\tt [cgxu,cgyu] = pdecgrad(p,t,c,U)} : sur le triangle num\'{e}ro $j,$
on a 
\begin{equation*}
\sigma (u)(x)=\left[ 
\begin{array}{cc}
\mathrm{cgxu}(1,j) & \mathrm{cgyu}(1,j) \\ 
\mathrm{cgxu}(2,j) & \mathrm{cgyu}(2,j)
\end{array}
\right] .
\end{equation*}
Et maintenant, \`{a} vous de jouer ! Vous rendrez un rapport de 4 ou 5 pages au
plus avec le r\'esultat de vos exp\'erimentations.
\end{itemize}

\end{document}