|
|
|
During year 2004/2005, i followed the DEA **analyse numerique** from the
|
|
|
|
university of [Paris VI](http://www.upmc.fr). This is now called
|
|
|
|
[parcours Analyse Numérique & Équations aux Dérivées
|
|
|
|
Partielles](http://www.ann.jussieu.fr/ANEDP).
|
|
|
|
|
|
|
|
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](http://en.wikipedia.org/wiki/Inpainting) in order to
|
|
|
|
predict the image at a lower lower. One application could be find new
|
|
|
|
compression algorithm. My advisors were [Albert
|
|
|
|
Cohen](http://www.ann.jussieu.fr/~cohen/) and [Simon
|
|
|
|
Masnou](http://www.ann.jussieu.fr/~masnou/).
|
|
|
|
|
|
|
|
[Read my report
|
|
|
|
(English)](http://www.freehackers.org/media/orzel/dea/stage/rapport.pdf)
|
|
|
|
|
|
|
|
## Videos
|
|
|
|
|
|
|
|
I did several videos that could have an academic interest.
|
|
|
|
|
|
|
|
[1](http://www.freehackers.org/media/orzel/dea/stage/Lenna_haar_reconstruction.avi)
|
|
|
|
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.
|
|
|
|
|
|
|
|
- [CompressionInterp-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-.5-.5.avi)
|
|
|
|
- [CompressionInterp-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-.5-1.avi)
|
|
|
|
- [CompressionInterp-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-.5.avi)
|
|
|
|
- [CompressionInterp-1-1.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-1.5.avi)
|
|
|
|
- [CompressionInterp-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-1.avi)
|
|
|
|
- [CompressionInterp-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-2.avi)
|
|
|
|
- [CompressionInterp-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1.5-1.avi)
|
|
|
|
- [CompressionInterp-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-2-1.avi)
|
|
|
|
- [CompressionInterp-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-2-2.avi)
|
|
|
|
- [CompressionLaplace-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-.5-.5.avi)
|
|
|
|
- [CompressionLaplace-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-.5-1.avi)
|
|
|
|
- [CompressionLaplace-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-1.avi)
|
|
|
|
- [CompressionLaplace-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-.5.avi)
|
|
|
|
- [CompressionLaplace-1-1.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-1.5.avi)
|
|
|
|
- [CompressionLaplace-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-2.avi)
|
|
|
|
- [CompressionLaplace-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1.5-1.avi)
|
|
|
|
- [CompressionLaplace-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-2-1.avi)
|
|
|
|
- [CompressionLaplace-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-2-2.avi)
|
|
|
|
- [CompressionDisoc-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-2-1.avi)
|
|
|
|
- [CompressionDisoc-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-2-2.avi)
|
|
|
|
- [CompressionDisoc-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1.5-1.avi)
|
|
|
|
- [CompressionDisoc-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-2.avi)
|
|
|
|
- [CompressionDisoc-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-.5-1.avi)
|
|
|
|
- [CompressionDisoc-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-1.avi)
|
|
|
|
- [CompressionDisoc-1-1.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-1.5.avi)
|
|
|
|
- [CompressionDisoc-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-.5.avi)
|
|
|
|
- [CompressionDisoc-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-.5-.5.avi)
|
|
|
|
|
|
|
|
## Horror story
|
|
|
|
|
|
|
|
On of the most frightening thing I've ever done with the computer. They
|
|
|
|
*requested* me to write something in
|
|
|
|
[Fortran](http://en.wikipedia.org/wiki/Fortran), so I did an
|
|
|
|
[implementation of the
|
|
|
|
heapsort](http://www.freehackers.org/media/orzel/dea/soft/tri_tas.f).
|
|
|
|
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
|
|
|
|
|
|
|
|
 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.`
|
|
|
|
|
|
|
|
- [Rapport du
|
|
|
|
projet](http://www.freehackers.org/media/orzel/dea/finance/rapport_finance.pdf)
|
|
|
|
- [Code source du projet
|
|
|
|
(bz2)](http://www.freehackers.org/media/orzel/dea/finance/sources-finance.tar.bz2)
|
|
|
|
- [Code source du projet
|
|
|
|
(tar.gz)](http://www.freehackers.org/media/orzel/dea/finance/sources-finance.tar.gz)
|
|
|
|
|
|
|
|
## 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
|
|
|
|
|
|
|
|
<!-- -->
|
|
|
|
|
|
|
|
- [Rapport du
|
|
|
|
projet](http://www.freehackers.org/media/orzel/dea/maillages/rapport.pdf)
|
|
|
|
- [Code source du projet
|
|
|
|
(bz2)](http://www.freehackers.org/media/orzel/dea/maillages/projetmaillage.tar.bz2)
|
|
|
|
- [Code source du projet
|
|
|
|
(tgz)](http://www.freehackers.org/media/orzel/dea/maillages/projetmaillage.tgz)
|
|
|
|
|
|
|
|
### 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é.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
### Projet particules
|
|
|
|
|
|
|
|
- [9 particules n'interagissant
|
|
|
|
pas](http://www.freehackers.org/media/orzel/dea/maillages/BouncingMolecule.avi),
|
|
|
|
elles rebondissent juste.
|
|
|
|
- [validation du tore pour les positions et
|
|
|
|
interactions](http://www.freehackers.org/media/orzel/dea/maillages/ValidationTore.avi)
|
|
|
|
(40 particules)
|
|
|
|
- [un phenomène intéressant ou des particules s'aggrègent puis
|
|
|
|
explosent
|
|
|
|
finalement.](http://www.freehackers.org/media/orzel/dea/maillages/agg.avi)
|
|
|
|
(40 particules)
|
|
|
|
- [une version ou tout
|
|
|
|
explose](http://www.freehackers.org/media/orzel/dea/maillages/comportementbizarre.avi),
|
|
|
|
pour montrer le type de bug/difficultés que l'on peut rencontrer
|
|
|
|
pendant le développement.
|
|
|
|
|
|
|
|
## Projet éléments Finis
|
|
|
|
|
|
|
|
(Cours de M.Perronnet)
|
|
|
|
|
|
|
|
L'ensemble du code du projet peut être télécharé ici:
|
|
|
|
|
|
|
|
- [Le projet
|
|
|
|
2D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem2d.tgz)
|
|
|
|
- [Le projet
|
|
|
|
3D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem3d.tgz)
|
|
|
|
|
|
|
|
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).
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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).
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Documents:
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffFirstTests.html)
|
|
|
|
- [Le fichier
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/orzel2d_2trous.cpp)
|
|
|
|
- [Le fichier mesh
|
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/orzel_section_2trous.msh)
|
|
|
|
|
|
|
|
### 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:
|
|
|
|
|
|
|
|
- [les modifications apportées a
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffTenseurConductivite.html)
|
|
|
|
- [Le fichier
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceConductivite.cpp)
|
|
|
|
|
|
|
|
### 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 :
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
*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:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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:
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffNeumann.html)
|
|
|
|
- [Le fichier
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceNeumann.cpp)
|
|
|
|
- [Le maillage grossier
|
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere.msh)
|
|
|
|
- [Le maillage fin
|
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere_big.msh)
|
|
|
|
|
|
|
|
### 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:
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffRobin.html)
|
|
|
|
- [Le fichier
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceRobin.cpp)
|
|
|
|
- [Le maillage fin
|
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere_big.msh)
|
|
|
|
|
|
|
|
### 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 :
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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*
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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:
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffProjet2d.html)
|
|
|
|
- [Le fichier
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceProjet2d.cpp)
|
|
|
|
- [Le maillage de la
|
|
|
|
pièce](http://www.freehackers.org/media/orzel/dea/piece.msh)
|
|
|
|
- [Le maillage de la pièce avec le
|
|
|
|
radiateur](http://www.freehackers.org/media/orzel/dea/piecerobin.msh)
|
|
|
|
|
|
|
|
### 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).
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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é... |
|
|
|
\ No newline at end of file |
