Capricelli/DEA

From Freehackers
Jump to: navigation, search

During year 2004/2005, i followed the DEA analyse numerique from the university of Paris VI. This is now called parcours Analyse Numérique & Équations aux Dérivées Partielles.


The last chapters are in French.

Research Internship

I worked on a new way of doing multiscale image analysis. The idea was to use the knowledege about disocclusion or inpainting in order to predict the image at a lower lower. One application could be find new compression algorithm. My advisors were Albert Cohen and Simon Masnou.

Read my report (English)

Videos

I did several videos that could have an academic interest.

[1] this one shows the image reconstruction from its decomposition on a Haar Basis. Haar basis are the first and simpler example of wavelets.

The following videos are experiments related to the methods described in my report.The two numbers are the corresponding values for alpha and beta.

Horror story

On of the most frightening thing I've ever done with the computer. They requested me to write something in Fortran, so I did an implementation of the heapsort. This was academic enought for them, and I'm now free to never use this programming language anymore.

Projet EDP et méthodes mathématiques en finance

Exempleput.png

Projet realisé dans le cadre de l'évaluation du cours de Messieurs Pironneau et Berestycki.

J'ai choisi le projet 10, sur les méthodes adaptatives pour la résolution de l'équation de Black-Scholes L'énoncé exact est :

 Run program 5.1 to 5.5. Adapt it to a barrier option with
 realistic data. If times allow make the change S -> S/(S+1)
 which maps R+ into (0,1) and compute a European put in this
 formulation and compare.

Projet sur les maillages adaptatifs

Projets realises dans le cadre de l'evaluation du cours de Messieurs Hecht et Frey, et réalisé en binôme avec Hanen Amor.

Nous avons choisit de faire deux projets:

  • Le projet 4 : Interpolation de maillages en O(n)
  • Le projet 8 : Schéma numérique et visualisation 3D pour un problème d'interaction de particules

Projet interpolation

Le test suivant vérifie que les lignes de courant sont conservées par notre programme d'interpolation. La première image montre une ligne de courant sur le maillage d'entrée. La seconde image montre la même ligne sur notre maillage interpolé.

Vit1 3d.001.png Vit0 3d.001.png

Projet particules


Projet éléments Finis

(Cours de M.Perronnet)


L'ensemble du code du projet peut être télécharé ici:

Le premier contient une copie de fem-v5, légérement modifiée comme indiqué dans les pages précédentes, les différentes versions du programme Laplace, et les fichiers de données correspondants.

Le deuxième contient l'ensemble du projet dans sa version 3D.


Premiers essais

J'ai modélisé la coupe transversale d'un tuyau transportant deux fluides à température constante. 0 pour le tuyau de gauche, et +10 pour le tuyau de droite. L'extérieur du tuyau est à température également constante, égale à +5. (On ne peut mettre que des conditions de Dirichlet pour l'instant).

Orzel2d 2trous maillage.png


Afin de mieux comprendre ce qu'il se passe, j'ai cassé la symétrie en imposant une source asymétrique dans l'épaisseur du tuyau. Pour ce premier exemple, le gradient est nul sur l'axe y=x, et est linéaire selon l'axe y=-x, négatif en bas à droite, positif en haut à gauche. (La formule est du type -x+y)

Orzel2d 2trous.png


La température sur l'ensemble de l'élément est majorée par 0 et 10, c'est-à-dire les températures extrêmes.

J'ai ensuite considéré un gradient nul sur l'axe y=-x, et linéaire selon l'axe y=x, négatif en bas à gauche, positif en haut à droite. (La formule est essentiellement x+y)

Contrairement à l'exemple précédent, la température dépasse les valeurs extrêmes, en négatif (bas-gauche) ou en positif (haut-droit).

Orzel2d 2trous2.png

Documents:

Prise en compte du tenseur de conductivité

Notre première modification vise à prendre en compte le tenseur de conductivité en x et y, c'est-à-dire en fonction du matériau, ou encore, en terme Mefisto, du Label.

On s'est contenté ici d'utiliser une matrice de conductivité diagonale, la multiplication par un vecteur revient donc à multiplier chaque composante par la composante correspondante de la diagonale.

L'équation différentielle vérifiée est alors -div[A(x)*grad(t(x))] = f(x), où t(x) représente la température en x, et A(x) la matrice du tenseur de conductivité. Dans la formulation variationnelle discrétisée, on remplace alors le calcul de l'intégrale sur un triangle élémentaire T de (grad(wi),grad(wj)) (produit scalaire), par le calcul de l'intégrale de (A(x)*wi,wj). Les wi représentent les fonctions chapeaux sur les noeuds du maillage.

Implémentation : on définit un tableau de matrices, l'indice du tableau sera le label du matériau, et chaque matrice sera en fait un tableau de deux éléments.

Documents:

Conditions de type Neumann

Jusqu'à alors, nous ne pouvions utiliser que des conditions de type Dirichlet. Je vais rajouter, dans un premier temps, la prise en compte de conditions de type Neumann. C'est le plus facile, car il ne faut modifier que le second membre.

Implémentation : Je parcours l'ensemble des côtés qui sont sur la frontière. Pour chaque côté, je vérifie d'abord qu'il fait partie d'une frontière où il y a des conditions de type Neumann. Sinon, bien sûr, je passe au côté suivant.

J'utilise alors une formule d'intégration numérique en 1D pour calculer l'intégrale sur le segment. Il ne reste plus qu'à modifier le second membre en conséquence.

Premier Test

Pour tester cette partie j'ai considéré ce premier exemple : Un carré avec un trou au milieu. La température est imposée (Dirichlet) sur le carré, avec un gradient linéaire selon l'axe y : -20 degrés en bas, et +20 en haut. J'ai imposé une condition de Neumann nulle sur le cercle intérieur, je dois donc obtenir des isothermes perpendiculaires à la frontière. Ce qui est bien le cas :

Condition neumann test.png

Tests

Je fais ensuite un test plus rigoureux. Je pars d'une solution connue, sur un objet facile à manipuler. Je vérifie que le résultat obtenu est bien celui attendu. J'ai pris comme objet un anneau, de cercle intérieur de rayon 4, et de cercle extérieur de rayon 10. La température est téta(x,y)=x*x+y*y. Ainsi, l'inverse du laplacien est constant égale à -4. Je considère alors le problème avec conditions de Dirichlet sur le cercle intérieur (température = 16) et condition de Neumann sur le cercle extérieur (le gradient vaut 2(x,y) et la normale unitaire au cercle vaut (x,y)/R2, donc la derivée normale vaut 2*R2). J'obtiens ceci:

Condition neumann test2.png

Ce qui est déjà prometteur, car la température est distribuée radialement entre 4*4 et 10*10.

J'ajoute un test en fin de programme qui calcule le max des différences entre les températures obtenues, et celles que j'aurais dû obtenir. Voici le code

 // test de validite
 R max = 0.;
 R tetamax =0.;
 assert (  Ug.N() == Th.nv );
 for ( int i = 0 ; i< Th.nv; i++ )
 {
         R diff;
         R2 P;
 
         P = (R2)(Th.vertices[i]);
         diff = Ug[i]-TetaExact(P) ;
         diff *= diff; // on prend le carre
         if ( diff>max ) max= diff;
         if ( Ug[i]>tetamax ) tetamax= Ug[i];
 }
 printf("Max des carres obtenu : %f\n", max);
 printf("Max des tetas obtenu : %f\n\n", tetamax);

Avec le premier maillage, correspondant à respectivement 16 et 32 segments sur les cercles intérieur et extérieur, j'obtiens :

 Max des carres obtenu : 0.750861  

Avec un maillage plus fin, (50 et 200 segments sur les cercles), j'obtiens :

 Max des carres obtenu : 0.005458  

Ce qui me parait raisonnable pour des températures entre 16 et 100.

Documents:

Conditions de type Robin/Fourier

L'étape suivante consiste à prendre en compte les conditions de type Robin/Fourier. Pour cela, il faut des données supplémentaires :

  • La température extérieure (fluide léchant le bord)
  • Le coefficient d'échange

Dans la formulation variationnelle correspondante, il faut alors modifier à la fois le second membre, et la matrice globale.

Pour le second membre, on fait exactement comme pour les conditions de Neumann, en remplaçant UNeumann(P) par GRobin(P) * TetaExtRobin(P) avec des notations évidentes.

Concernant les modifications portant sur la matrice, il faut intégrer g.wi.wj, pour chaque segment sur la frontière de Robin. Si je considère g constant, je peux utiliser la formule intégrale(f) sur [0,1] ~ 1/6(f(0)+4f(1/2)+f(1)), qui est exacte sur les polynômes P2 que je considère ici. Après changement de variable pour se ramener au segment courant, cela revient à g*longueur(segment)/6, si i est différent de j, et à g*longueur(segment)/3 si i==j.

Modifications apportées à fem v5

Pour modifier la matrice globale, je n'ai pas trouvé de fonction correspondante dans l'API de fem, j'ai donc rajouté la fonction Add2x2Matrice aux fichiers MatriceCreuse_tpl.hpp et MatriceCreuse.hpp, dont voici le code (clairement inspiré de la fonction operator+=)

 template<class R>
 MatriceCreuse<R> & MatriceProfile<R>::add2x2Matrice
         (int i0, int i1,R a00, R a01, R a11) {
   int i,j;
   if (!D) { // matrice vide·
     cerr << "Big bug empty matrice not handled (yet) in add2x2Matrice" << endl;
     exit(1);
   }
   D[i0] += a00;
   D[i1] += a11;
 
   assert( i0 != i1 );
 
   i=i0; j=i1;
         if      (j < i)  L[ pL[i+1] - (i-j) ] += a01;
         else if (j > i)  U[ pU[j+1] - (j-i) ] += a01;
 
   i=i1; j=i0;
         if      (j < i)  L[ pL[i+1] - (i-j) ] += a01;
         else if (j > i)  U[ pU[j+1] - (j-i) ] += a01;
 
   return *this;
 }

Tests

Pour tester cette partie, j'ai utilisé le même exemple qu'avec les conditions de Neumann. J'impose toujours la température sur la frontière intérieure (Dirichlet, donc). Je fixe une température extérieure (TetaR=1000 dans cet exemple). La valeur de g se calcule facilement : g = - 2*R2/(TetaR-R2*R2).

J'obtiens, avec le maillage fin :

 Max des carres obtenu : 0.299266

Qui paraît, encore une fois, raisonnable. Je note quand même que la précision est moindre que dans l'exemple avec Dirichlet+Neumann.

Documents:

Exemple concret en 2D, sans condition de robin

Pour traiter un exemple, j'ai choisi de modéliser la température d'une pièce en L dont voici le maillage :

Exemple2d sans robin maillage.png

L'intérieur représente l'air de la pièce, les trois petites surfaces représentent des murs (oui, assez épais, disons qu'il s'agit d'un morceau de château..). Le bâtiment se prolonge sur la gauche et vers le bas, c'est pour cela que je n'ai pas mis de murs à ces endroits. Il y a deux fenêtres, une grande et une petite. Je n'ai pas modelisé de portes car c'est inutile vu les conditions utilisées.

Bien entendu, l'air de la pièce et les murs n'ont pas la même conductivité.

Les conditions imposées sont les suivantes :

  • Température de 0 degré à l'extérieur, sur le mur et la fenêtre de droite (Dirichlet)
  • Température de 20 degrés sur les murs intérieurs (Dirichlet)
  • Le mur et la fenêtre du haut sont baignés par le soleil (Condition de Neumann)

Voici ce que j'obtiens:

Exemple2d sans robin result.png

On note que le gradient est bien plus grand dans les murs, qui jouent évidemment le rôle d'isolant par rapport à l'extérieur. On observe bien de quelle manière la fenêtre de droite fait passer le froid.

Exemple complet avec condition de robin

Exemple2d avec robin maillage.png

Exemple2d avec robin result.png

La température est alors plus homogène dans la pièce (Regarder la surface où la température est entre 20 et 18 degrés)

Documents:


Passage à la 3D

Il s'agit maintenant de faire la même chose, mais en 3D avec des tétraèdrisation P1-Lagrange. C'est loing d'être une petite modification locale à un fichier et un minimum de méthodologie était nécessaire.

Modifications portant sur le projet

Dans un premier temps, j'ai recopié l'ensemble du répertoire, et effacé tous les fichiers dont je n'avais pas besoin (comme par exemple Element_P2h.cpp ou BamgFreeFem.*). Le fichier Makefile est modifié en conséquence.

Puis j'ai créé les classes R3 (déjà à peu près faite) et Tetraedre. J'ai modifié l'existant pour qu'au lieu d'utiliser des points R2, des segments et des triangles, on utilise respectivement des points R3, des triangles et des tetraèdres. En particulier, j'ai trouvé pratique de ne pas du tout utiliser R2, et de refaire R3 from scratch. La première raison est que R3 héritait de R2 et cela soulevait de nombreux bugs (par exemple si l'opérateur produit scalaire n'est pas défini dans R3, le compilateur caste automatiquement et silencieusement les points sur des R2 et utilise le produit scalaire R2..). La deuxième raison est que cela permet, grace aux erreurs du compilateur, de trouver toutes les occurences de R2 à remplacer.

Ensuite j'ai modifié la fonction Mesh::read() pour être capable de lire les fichiers générés par Mefisto dans le cas 3D. J'ai utilisé le fichier hexaedre.msh pour valider cette partie.

S'est alors posé le problème de FESpace. Cette classe, ainsi que celles associées, sont difficiles à comprendre et surtout à modifier. Ainsi, après plusieurs jours passés à me battre avec ces fichiers, et ayant constaté que ce dont je me servais dans FESpace était en fait directement issu de la classe Mesh, j'ai décidé de supprimer cet intermédiaire. Ce qui n'a posé aucun problème (ce qui ne signifie pas pour autant que cela ait été rapide). Par exemple le champs FESpace.NbOfDF est en fait TTh.nv, etc...

Attention, la classe FESpace a un rôle dans fem-v5 En particulier, mais pas seulement, celui de gérer les Mortars. C'est parce que je limitais les possibilités du programme que j'ai pu faire cela.

Pour la fonction gradgrad(), j'ai besoin du gradient des fonctions chapeau sur les tétraèdres. Si on considère un sommet A d'un tétraèdre, alors, le gradient de la fonction chapeau valant 1 en A et 0 sur les autres sommets est un vecteur, constant sur le tetraedre, de direction la hauteur issue de la face opposée a A, et de longueur l'inverse de cette hauteur. Ce qui est realisé par la fonction

 R3 Tetraedre::H( int i )  const
 {  assert( i>=0 ); assert( i<4 );
 
   // si i est impair, on inverse c et b pour
   // avoir les indices dans l ordre direct
   R3     &D = *(vertices[  i      ]),
          &C = *(vertices[ (i+1 + (i%2))%4 ]),
          &B = *(vertices[ (i+2 - (i%2))%4 ]),
          &A = *(vertices[ (i+3)%4 ]);
 
   // check we took the coordinates in the right order
   assert  ( ((B-A)^(C-A),(D-A)) >0 );
 
   // On cherche HD, avec H dans le plan ABC et HD perpendiculaire a ce plan
   R3 h = ( B-A )^( C-A ); // normal a la face ABC
   R norm = Norme2(h);;
   h /= norm; // h est unitaire
   h /= (6.* volume/norm); // on divise par la longueur de la hauteur
   return h;
 } 

J'ai reimplémenté les formules d'intégrations numériques sans utiliser celles de fem v5. Par curiosité, et aussi dans un soucis de clarté. En effet, fem-v5 utilise trois fichiers (FormuleIntegration.cpp FormuleIntegration.hpp QuadratureFormular.hpp), les coefficients sont normalisés (QuadratureFormular1d) ou pas (P_QuadratureFormular_T*), l'implémentation utilise parfois des structures, parfois des tableaux ou même des classes..

J'ai réussi à garder le fichier gibbs.cpp, qui contient un algorithme pour renumeroter les éléments du maillage pour optimiser le skyline (et donc la taille en mémoire) de la matrice creuse.

Modifications portant sur le programme principal

Quelques modifications triviales : les matrices diagonales des tenseurs de conductivités ont maintenant 3 éléments, passage de area à volume dans la fonction gradgrad(), on utilise plus FESpace..

Puis, pour chaque étape, on adapte à la 3D, en se basant sur notre nouvelle infrastructure. Pour construire le second membre, on ne travaille plus avec des triangles mais avec des tétraèdres et on utilise une formule d'intégration 3D. Et Ainsi de suite...

Tests

Vu l'ampleur des modifications apportées, j'ai prévu de nombreux tests, du plus simple au plus complet. Pour chaque test, j'ai un fichier (par exemple donnees_dirichlet_uniforme.cpp), et selon le test que je veux effectuer, j'inclus le fichier correspondant. Cela evite les multiples #if/#endif dans le code et le rends plus lisible.

Le premier test consiste à imposer une température nulle sur toutes les faces. On test ainsi que le programme compile et link, et qu'au moins sur un exemple aussi trivial, on obtient bien la bonne température.

Le deuxieme test (donnees_dirichletnul.cpp) impose une température nulle au bord, mais avec un fOmega non nul. Il permet de tester la construction du second membre, la construction et l'inversion de la matrice.

Le troisième test est une fonction quelconque. J'ai pris 2.x.y.z-10.y.x^2+6.x.z^2 , qui n'est pas nulle sur le bord, et dont le laplacien n'est pas nul non plus : cela permet de tester, en plus des tests précédents, les conditions au bords : Dirichlet, puis Neumann, et enfin Robin/Fourier.

Enfin le test donnees_x2y2z2.cpp est l'exemple classique, et m'a permis de comparer les resultats obtenus avec d'autres eleves.

Résultats

  • donnees_dirichlet_uniforme.cpp : j'obtiens exactement la temperature attendue.
  • donnees_dirichletnul.cpp : la aussi j'obtiens exactement la temperature attendue.
 Max des carres obtenu : 0.000000
 sqrt   : 0.000000
 Min,Max des tetas obtenu : 0.000000, 1.926837
 Min,Max des tetas reels  : 0.000000, 1.926837
 Taux d'erreur : 0.000000
  • donnees_quelconque.cpp : donnees de Dirichlet : La aussi j'obtiens exactement la temperature attendue.
 Max des carres obtenu : 0.000000
 sqrt   : 0.000000
 Min,Max des tetas obtenu : -100.000000, 60.000000
 Min,Max des tetas reels  : -100.000000, 60.000000
 Taux d'erreur : 0.000000
  • donnees_quelconque.cpp : donnees de Dirichlet, de Neumann et de Robin/Fourier : (c'est le cas le plus difficile)
 Max des carres obtenu : 62.684211
 sqrt   : 7.917336
 Min,Max des tetas obtenu : -96.800396, 53.333334
 Min,Max des tetas reels  : -100.000000, 60.000000
 Stats under,toomuch   : 830, 54
 Taux d'erreur : 4.948335


Exemple 3D

Pour l'exemple 3D, j'ai voulu modéliser la température de l'eau d'un chauffe-eau, tel que ceux utilisés pour faire du thé, du café.. J'ai donc modélisé sous Mefisto un cylindre (le contenant), et un hexaèdre oblong (résistance chauffante).

Exemple3d.cylindre.png

Exemple3d.resistance.png

J'ai triangulé ces deux surfaces..., j'ai 'amelioré' (option 30 du menu surfaces) ces deux surfaces. J'ai ensuite créé un volume (option 8) avec ces deux surfaces fermées. Puis j'ai essayé de tétraedriser ce volume grace à l'option 9 du menu surface. Je n'ai jamais reussi. Je me suis battu pendant une semaine, et jusqu'a 5h du matin la veille de la soutenance mais je n'obtient que "TABLEAU PTXYZD SATURE".

J'ai essayé avec l'hexaèdre seul, j'ai essayé avec une sphère, j'ai essayé tous les paramètres possibles pour cette option 9. C'est un echec.

UPDATE M.Perronnet m'a expliqué le problème : la surface inférieur de l'hexaèdre de la résistance est confondu avec le disque inférieur du cylindre, c'est ce qui pose problème à Mefisto. A la rigueur, si le maillage correspondait, cela serait peut-être passé...