\documentclass[11pt]{article}
\usepackage{graphicx, algorithm,pdfsync, bm, amsfonts, amsmath,color}


\usepackage[latin1]{inputenc}
\usepackage[french]{babel}
\usepackage{hyperref}
\usepackage{FFF}
\def\textRed{\color{red}}
\def\textBlack{\color{black}}
\newcounter{thequestion}
\setcounter{thequestion}{0}
\def\Question#1{ \addtocounter{thequestion}{1}\par\noindent {\bf Question \arabic{thequestion} {#1}} }
\def\p{\partial}
\def\d{\hbox{d}}
\def\R{\mathbb{R}}
\title{Projet Informatique Scientique (MM031)  2011-12: Surface Minimale, v3 }

\author{Fr\'ed\'eric Hecht}

\def\d{\partial}
\begin{document}
\maketitle

\textbf{Projet2 : Il est possible de faire le projet en binôme,  il  est à rendre  le vendredi 11 mai 2012  24h00 par mail à Frederic.hecht@upmc.fr, 
 (faire un  archive  . tar.gz, ou .zip  des sources et rapport) }

\textbf{Projet1  Attention les modification  par rapport à  la version original sont en {\color{red} rouge}, et de plus il est possible de faire le projet en binôme
qui est à rendre  le 28 février 2012  24h00 par mail à Frederic.hecht@upmc.fr }



\section{G\'{e}n\'{e}ralit\'{e}s}
\paragraph{G\'eom\'etries}
On s'int\'{e}resse au calcul du problème de surface minimale et à sa visualisation.

Le  problème de surface minimale  est:  soit $\Gamma$ une courbe fermée de $\R^3$, le but est de trouver la surface ${S}$ d'aire 
minimale de bord  la courbe $\Gamma$.

Premièrement, nous supposerons que la surface est paramètré par $F_u : (x,y) \mapsto (x,y,u(x,y)) $ et que $\Gamma$  sera l'image
du bord $\d\Omega$ d'un domaine  $\Omega$ de $\mathbb{R}^2$ par $F_u$, c'est à dire que $u$ est connue sur $\d\Omega$ et l'on
écrira que $ u = g $ sur $\d\Omega$ où $g$ est une fonction connue continue.

Maintenant notre problème peut s'écrire :
Trouver la fonction $u$ tel que 
\begin{equation}
 u = \mathop{\mathtt{arg\ min}}_{u_{|\p\Omega}=g} \; J(u)  
\end{equation}
où
\begin{equation}
J(u)= \int_{\Omega} \sqrt{1+|\nabla u|^2} d\Omega
\end{equation}

\section{Projet 1}

Dans un premier temps nous voulons étudiez ce problème en dimension 1,
C'est a dire que l'on va supposer que u est invariant
 par toutes rotation de centre l'origine
 


%\Question{}


Nous nous proposons d'{\'e}crire une m{\'e}thode num{\'e}rique pour r{\'e}soudre le probl{\`e}me non-lin{\'e}aire suivant :

\noindent Trouver $u$ une fonction de $]a,b[\mapsto \mathbb{R}$, avec $0<a<b$, 
qui v{\'e}rifie les conditions aux limites  (CL)
$$ u(a)=\alpha, \qquad u(b)=\beta, \quad \mbox{(avec $\alpha$ et $\beta$ des valeurs donn{\'e}es)},$$
et qui r{\'e}alise {\bf le minimum} de la fonctionnelle $J$ suivante :
$$ J(u) = {2\pi} \int_a^b r\, \sqrt{1+ (u'(r))^2} \,dr. \eqno (P)  $$

\noindent On supposera que ce probl{\`e}me est bien pos{\'e} dans $H^1(]a,b[)$, c'est {\`a} dire, il existe une unique  solution $u$
et qui ne  d{\'e}pend contin{\^u}ment des 4 donn{\'e}es $a, b, \alpha, \beta$.\\

\begin{description}
\item[Q0] Montrer que $J$ est la surface (de r{\'e}volution) d{\'e}finie (balay{\'e}e) par la rotation de la courbe $z=u(r)$, pour $ a < r < b $, autour de l'axe  $z$, o{\`u} $(r, z)$ sont les coordonn{\'e}es  polaires, avec $r=\sqrt{x^2+y^2}$.\\

\item[Q1] Calculer la diff{\'e}rentielle de Fr{\'e}chet  $DJ(u)$ de $J$ en $u$ ; $DJ(u)$ est une application lin{\'e}aire de $H^1_0(]a,b[)$ dans $\mathbb{R}$.\\



\item[Q2] Justifier pourquoi $u$ est solution du probl{\`e}me suivant :
$$   \forall v \in H^1_0((]a,b[), \quad {2\pi} \int_a^b r \frac{u'v'}{ \sqrt{1+ (u'(r))^2} }dr = 0$$\\

\item[Q3] Montrer que la cat{\'e}no{\"\i}de  $u(r) = acosh(r)$ (i.e. $cosh(u)=r$) est la solution du probl{\`e}me pr{\'e}c{\'e}dent pour $a=1$, $\alpha=0$, et {\color{red} $ \beta= acosh(b)$ }. \\


\item[Q4] On veut discr{\'e}tiser le probl{\`e}me (P) avec des {\'e}l{\'e}ments finis $P_1$ Lagrange : soit $N$ un entier strictement positif,
on d{\'e}finit $h=(b-a)/N$ , $q_i = a + hi $ , pour $i = 0, .., N$, les points du  maillage et l'espace $V_N$ des fonctions de $\mathcal{C}^0(]a,b[)$  affines par morceaux sur $]q_i,q_{i+1}[$ pour $i=0, .., N-1$.

Soit $w_i$ la base de $V_N$ telle que $ w_i(q_j) = \delta_{ij}$, o{\`u} $\delta_{ij}$ est le symbole de Kroneker.

Soit une fonction $v_N$ d{\'e}finie par ses valeurs aux n{\oe}uds du maillage $$v_N= \sum_{i=0}^{N} v_i w_i,\quad  \mbox{avec} \quad   (v_i)_{i=0}^N \in \mathbb{R}^{N+1}.$$
Ecrire un algorithme formel qui calcule  la fonctionnelle $J(v_N)$
en fonction des valeurs $v_i$.\\

%\item[Q5] Ecrire une base de l'espace vectoriel $V_N\cap H^1_0(]a,b[)$.\\

\item[Q6N] On veut utiliser l'algorithme de minimisation {\`a} pas fixe  qui s'{\'e}crit comme suit :\\
Soit $ u^0_N$ un {\'e}l{\'e}ment de $V_N$ donn{\'e}, v{\'e}rifiant les conditions aux limites (CL) ;\\  on calcule r{\'e}cursivement
$$ u^{n+1}_N = u^{n}_N - \rho \nabla J(u^{n}_N),$$ o{\`u} $\rho$ est une constante d{\'e}finie par l'utilisateur, et le gradient $ \nabla J(u)$ est d{\'e}fini  {\`a} partir d'un produit  scalaire $(,)_N$ de $V_N$ par :
$$ (\nabla J(u),v)_N = DJ(u)(v). $$ .  On prendra comme produit scalaire  

\begin{enumerate}
\item   Soit  $$(u,v)_N = \sum_{i=0}^N u_i v_i,\quad \mbox{o{\`u}}\quad
v= \sum_{i=0}^N v_i w_i,\quad \mbox{et}\quad u= \sum_{i=0}^N u_i w_i.$$
\item  ou soit $$(u,v)_H = \int_\Omega u' v '.$$
\end{enumerate}
Evaluer les composantes $ g^n_i$ du gradient  d{\'e}finies pour   le produit scalaire $(,)_N$  par   $$\nabla J(u^{n}_N) = \sum_{i=0}^N g^n_i w_i,$$ comparer les par rapport 
${\partial J(u^{n}_N)}/{\partial u_i}$, avec $ u^{n}_N= \sum_{i=0}^N u_i w_i.$  et  {\'e}crire compl{\`e}tement l'algorithme  composante par composante.

\item[Q6H] Comment  calculer le gradient  $\nabla J(u^{n}_N)$ d{\'e}finies pour le produit scalaire  $(,)_H$. 
\end{description}



%\section{Partie pratique C++  (1 heure)}

Le but de cette partie est de programmer l'algorithme pr{\'e}c{\'e}dent.

On supposera que cette d{\'e}finition est active :\verb! typedef double R; !


\begin{description}
\item[Q7]   Ecrire une fonction \texttt{C++}  qui
calcule l'aire d'un tronc de c{\^o}ne $c(r)$ , de prototype  \verb!AireTroncDeCone(R ri, R ri1, R vi, R vi1) ; !   d{\'e}fini pour   $ \mathrm{ ri} < r <\mathrm{ri1}$ ; $c(r)$ la fonction  affine  telle que $c(\mathrm{ri})=\mathrm{vi}$
 et $c(\mathrm{ri1})=\mathrm{vi1}$.\\

\item[Q8] Ecrire une fonction    qui {\'e}value $J$ pour la fonction $v= \sum_{i=0}^N v_i w_i$, repr{\'e}sent{\'e}e par le vecteur \texttt{v} de ses valeurs, et qui a  pour prototype de \texttt{C++} \\
 \verb!R J(R a, R b, int N, R *v) ; !.\\

\item[Q9] Ecrire une fonction   qui {\'e}value $\nabla J$ pour la fonction $v= \sum_{i=0}^N v_i w_i$
et qui retourne le pointeur $g$ qui est un argument qui va  stocker le vecteur $\nabla J$,
et qui a pour prototype \texttt{C++}   \verb!R * GJ(R a, R b,int N,R *v,R *g) ; !

\item[Q10N] Supposons que  l'algorithme de la question Q6 est programm{\'e} avec  la fonction suivante :

\bFF
 @bool AlgoQ6N(R a,R b,R alpha, R beta, VN & un,R eps,
              \textRed R rho \textBlack,int itermax=10000); // ajoute de rho 
 {
    @int N= un.N;
    VN gn(N),
    un(0) = alpha;
    un(N) = beta;
    @for(@int n=0;n<itermax;++n)
      {
         gn = GJ(a,b,N,un,gn);
         \textRed R err= gn.linfty();
         \textRed cout << " n = "<< n << " " << err<< endl; 
          @if(err  <eps) @return true;
         \textRed  gn *= rho;  //  changement ICI  
         un -= gn; 
      }
      @else @return @false;
 }
\eFF
Ecrire une classe tableau \texttt{VN} qui dérive de \verb!vector<double>! ainsi que les prototypes  des m{\'e}thodes et des fonctions associ{\'e}es qui seront n{\'e}cessaires  pour la fonction \verb! AlgoQ6! fonctionne. Tester et valider cette fonction sur l'exemple de la question Q3.

Afficher le solution calculé et la solution exact avec gnuplot. 

\item[Q10H] Même question que Q6N mais avec le gradient calculer avec le produit scalaire $(,)_H$

\item[Q10W] Ecrire la méthode de Newton pour résoudre le problème non linéaire suivant: 
Trouver $u \in V_N$ , tel que $u(a)=\alpha, u(b)=\beta$  tel que

$$   \forall v_N \in V_N\cap H^1_0(]a,b[) , \quad {2\pi} \int_a^b r \frac{u'v'}{ \sqrt{1+ (u'(r))^2} }dr = 0$$



\item[Q11] Comparer les résultats, et les temps de calcules entre les algorithmes des questions Q10N,Q10H,Q10W. 
 \item[Q12] Faire une étude de convergence des algorithme en fonction de $N=3,10^5$. 

%Rappel : les r{\`e}gles de conversion (\og cast\fg\ en Anglais)  d'un type $T$ en $A$  sont g{\'e}n{\'e}r{\'e}es, soit {\`a} partir d'un constructeur \texttt{A(T)}
%    dans la classe $A$, soit en utilisant l'op{\'e}rateur de conversion  \texttt{operator (A) () } dans la classe $T$ (ou
%     \texttt{operator (A) (T)} hors de la classe). \\
%     %De plus il ne faut pas oublier  que le compilateur \texttt{C++ } fait automatiquement au maximum  un niveau de conversion pour trouver la bonne fonction ou  le bon op{\'e}rateur.\\
%
%\item[Q11] Quel est le code (op{\'e}rateurs, constructeurs,...) qu'il faut ajouter dans la classe \texttt{VN}
%    pour {\'e}liminer la possibilit{\'e} de perdre de  m{\'e}moire (toutes les allocations faites avec un \texttt{new} ont un \texttt{delete} correspondant)?
%    %Bien sur il faut qu'il n'y ai aucun fuite de m{\'e}moire de cette classe, qu'elle est les code qu'il faut ajouter\\
%
%\item[Q12] Ecrire une programme qui teste et valide votre algorithme num{\'e}rique.
\end{description}




\section{Projet 2}
\def\Th{\mathcal{T}_h}
 
Ce problème n'est pas simple théoriquement, mais si on le regarde d'une point de vue discret, il relativement simple.

Nous allons supposer que nous disposons d'un maillage $\Th$ de notre domaine plan $\Omega$, et soit
$ V_h = \{ v \in H^1(\Omega) / \forall K \in \Th ; v_{|K} \in P1(K) \} $ 
soit $u_h \in V_h$, la surface  facetiser que l'on veut minimise est $\sigma_{u_h} = F_{u_h}(\Omega)$
où $ F_{u_h} : (x,y) \mapsto (x,y, u_h(x,y)$, pour que notre problème est un sens physique et ne dégénère pas trivialement,

$u_h$ sera donne égal à $\Pi_h(g) $ sur $\Omega$, où $\Pi_h$ est interpoler dans $V_h$ de la fonction $g$ donné. 

 


Donc notre problème  est de trouver 
 
$$ \mathop{\texttt{arg min}}_{u_h\in V_h, {u_h}_{|\p \Omega} =\Pi_h(g) } J_1=  { \sum_{K\in \Th} |F_{u_h}(K)| } \eqno (P_1)
$$
 
Maintenant, on peut prendre une paramétrisation plus générale du type $ U_h = (u_{1h}, u_{1h}, u_{1h}) \in V_h^3$ 
ou l'on fait bouge les $3$ composante du maillage

$$ \mathop{\texttt{arg min}}_{U_h\in V_h, {U_h}_{|\p \Omega} =\Pi_h(G) } J_3=  { \sum_{K\in \Th} |{U_h}(K)| } \eqno (P_3)
$$
 
les domaines de calcules seront $ \Omega_1$: une couronne entre $r_a=a$ et $r_b=b$, ou un rectangle  $\Omega_2 =]x_0,x_1[\times ]y_0,y_1[$.

De plus Parmi les condition aux limites nous  voulons étudions  le cas  correspondant 
où sur le bord du maillage la surface est défini par 
$ x,y, \lambda*sin(\pi 2)*cos(\pi x)$ où $\lambda$ est un paramètre. 


Partie Freefem++  
\begin{description}
\item[Q0]  Petite recherche bibliographique afin de pouvoir valider votre projet. 
\item[Q1] Montrer que $$J_1(u_h) = \int_{\Omega} \sqrt{1+ \nabla u.\nabla u } dx$$
et trouver $DJ_1(u_h) v_h $ , implementer en FreeFem++ la méthode de gradient à pas simple 
pour minimiser $J_1(\Omega)$, valider avec les projet 1, sur $Omega_1$ en choisissant la fonction $g$ et $a$ et $b$
pour des $ h = 0.2,0.1,0.05$. De plus il faut  choisir comme produit scalaire 1) le produit scalaire de de $\mathbb{R}^n$,
puis le produit scalaire de $H^1_0(\omega)$ , $(u,v) = \int_\Omega \nabla u.\nabla v $. 

\item[Q2] Montrer que $$J_2(u_h) = \int_{\Omega} \sqrt{||\p_x U \wedge  \p_y U || } dx$$
et trouver $DJ_2(u_h) v_h $ , implementer en FreeFem++ la méthode de descente à pas simple 
pour minimiser $J_1(\Omega)$, valider avec $U=(ycos(x),ysin(x),x)$ sur $\Omega_2=[0,4\pi]\times]0,4\pi[$ en choisissant la fonction .
pour une suite de trois maillages raisonnable.  

\item[Q3] Utiliser le méthode \texttt{NLGC} de freefem++ avec ou sans preconditionnement pour minimiser $J_1$

\item[Q4] Utiliser le méthode \texttt{NLGC} de freefem++  sans preconditionnement pour minimiser $J_3$
\item[Q5] Ecrire la méthode Newton pour Minimiser $J_1$ et $J_3$ en utilisant quel itération  des algorithme précendant pour
initialer la méthode Newton, si il y a des problème de  convergence. 

\item[Q6] calcule numériquement  en $\lambda = 1$ , $$\frac{\p \mathop{\texttt{arg min}}_{u_h\in V_h, {u_h}_{|\p \Omega} =\Pi_h(g) } J_1}{\p \lambda} \eqno (dJ1)$$  et 
$$\frac{\p\mathop{\texttt{arg min}}_{u_h\in V_h, {u_h}_{|\p \Omega} =\Pi_h(g) } J_3}{\p \lambda} \eqno (dJ3)$$ pour la 
surface de bord défini par $ (x,y, \lambda*sin(pi*2)*cos(pi*x)) $, sur un maillage du carré unité $20\times20$.




\item[Q7] Un utilisant les classes du Cours,
Implementer en \texttt{C++}  dans un programme 
qui est dynamiquement parametrable, via \texttt{argc} et \texttt{argv}, les paramètres du \texttt{shell}
passer a la fonction \texttt{main}, 
c'est-à-dire \verb! int main(int argc, const char** argv) !.

les paramètres seront
le maillage, puis l'algorithme : 
\begin{enumerate}
 \item   Algorithme de Gradient à pas fixe non précondionné pour $J_1$
 \item   Algorithme de Gradient à pas fixe  précondionné pour $J_1$
 \item   Algorithmes de Gradient à pas fixe non précondionné $J_3$
 \item   Algorithme de Gradient à pas fixe non précondionné pour $J_1$
 \item   Algorithme de Newton pour  $J_1$
 \item   Algorithme de Newton pour $J_3$
\end{enumerate}

plus des paramètres telle que:  choix de rho, nombre itération, test d'arrêt, choix des conditions aux limites.

\item[Q8] Ajouter le calculer  de $(dJ1)$ et $(dJ3)$ par différence fini et  par différentiation
automatique.
\item[Q9] Mettre ces algorithmes dans \texttt{GLplot.cpp} telle que les itérations soit faite avec
des timer GLUT, et où les paramètres sont changé dynamiquement dans des menus. 

\item[Q10] Pour finir, un petit interface java graphique  du programme 

\end{description} 
 
%
%
%
%\item[Q2]  En utilisant les timer   de Glut (via  \texttt{glutTimerFunc}) écrire un algorithme
%de descente à pas fixe pour minimiser $S_h$ avec la paramétrisation i) avec un gradient calculer par
%différence finite. Le timer fera une itération. 
%
%\item[Q3]  En utilisant les timer   de Glut (via  \texttt{glutTimerFunc}) écrire un algorithme
%de descente à pas fixe pour minimisez $S_h$ avec la paramétrisation ii) avec un gradient calculer par
%difference finite. Le timer fera un itération. 
%\item[Q4] Utiliser les outils de la différenciation automatique pour calculer le gradient de $J$ avec les deux paramétrisations. 
%
%\item[Q5] Utiliser un algorithme de type gradient conjugué non linéaire préconditionné par un laplacien avec la paramétrisation i).
%\item[Q6] Utiliser un algorithme de type gradient conjugué non linéaire préconditionné  avec la paramétrisation ii).
%\item[Q7] Comment validez vos algorithmes?
%\end{description}

%\bigskip
%A faire:
%
%\begin{enumerate}
%\item  Pour le jeudi 1 décembre 2011, envoyé votre binôme par mail \url{mailto:frederic.hecht@upmc.fr?subject=BINOME%20MPE} 
%\item  Le Dimanche 4 décembre, je met sur le web, liste de binôme et type de projet à faire. 
%\item Pour le mardi 13 décembre 2011, envoyé (1 par binôme) votre \texttt{class R3}, et la fonction 
%\begin{verbatim} 
%double AireDF(const R3  P[3], R3 DP[3] ) ;  
%\end{verbatim}
%qui retourne l'aire $s$ du triangle  de sommet P[0], P[1], P[2] et
%et où  \[\mathtt{DP[i]}.\xi  \sim   \frac{\partial s}{\mathtt{P[i]}.\xi} , \quad \xi \in \{ x,y,z\}, \; i \in \{0,1,2\}\]
%et la dérive partiel  est calculer via une méthode au différence. 
%\item Pour écriture vos classe vous  utilisez les idée qui sont définies dans \url{http://www.ann.jussieu.fr/~hecht/ftp/MPE/femgrid2d.tgz}
%\item En utilisant les outils de différentiation  automatique 
%recoder la fonction \texttt{AireDF} en: 
%\begin{verbatim} 
%double AireDFAD(const R3  P[3], R3 DP[3] ) ;  
%\end{verbatim}
%
%\item Finir remarquer que DJ et J peuvent être calcul avec les code suivant (non testé)
%\begin{verbatim} 
%  KNM<double> DJ(Th.nv,3); 
%  double J=0;
%  DJ=0. ; 
%  R3 P[3], DP[3];
%  for(int k=0;k<Th.nt;++k)
%  	{
%		Triangle K=Th[k];		
%	    P[0]=K[0];
%	    P[1]=K[1];
%	    P[2]=K[2];
%	    J+= AireDF(P,DP);
%	   
%	    for(int i=0;i<3;++i)
%	    {
%	    	  if(K[i].onGamma()) continue;
%		   DJ(Th(k,i),0) += DP[i].x;
%		   DJ(Th(k,i),1) += DP[i].y;
%		   DJ(Th(k,i),2) += DP[i].z;
%	    }  
%	} 
%\end{verbatim}
%
%\item Pour  le vendredi 6 Janvier, écrire une classe correspondant votre projet à partir des exemples 
%disponible dans \texttt{femgrid2d.tgz} et visualisez (Question Q1). 
%
%\end{enumerate}
%
%
%




  




\end{document}
