|
|
During year 2004/2005, i followed the DEA **analyse numerique** from the
|
|
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).
|
|
|
university of [Paris VI](http://www.upmc.fr). This is now called
|
|
|
|
|
[parcours Analyse Numérique & Équations aux Dérivées
|
|
The last chapters are in French.
|
|
|
Partielles](http://www.ann.jussieu.fr/ANEDP).
|
|
|
|
|
|
|
## Research Internship
|
|
|
The last chapters are in French.
|
|
|
|
|
|
|
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/).
|
|
|
## Research Internship
|
|
|
|
|
|
|
[Read my report (English)](http://www.freehackers.org/media/orzel/dea/stage/rapport.pdf)
|
|
|
I worked on a new way of doing multiscale image analysis. The idea was
|
|
|
|
|
to use the knowledege about disocclusion or
|
|
## Videos
|
|
|
[inpainting](http://en.wikipedia.org/wiki/Inpainting) in order to
|
|
|
|
|
predict the image at a lower lower. One application could be find new
|
|
I did several videos that could have an academic interest.
|
|
|
compression algorithm. My advisors were [Albert
|
|
|
|
|
Cohen](http://www.ann.jussieu.fr/~cohen/) and [Simon
|
|
[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.
|
|
|
Masnou](http://www.ann.jussieu.fr/~masnou/).
|
|
|
|
|
|
|
The following videos are experiments related to the methods described in my report.The two numbers are the corresponding values for alpha and beta.
|
|
|
[Read my report
|
|
|
|
|
(English)](http://www.freehackers.org/media/orzel/dea/stage/rapport.pdf)
|
|
- [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)
|
|
|
## Videos
|
|
- [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)
|
|
|
I did several videos that could have an academic interest.
|
|
- [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)
|
|
|
[1](http://www.freehackers.org/media/orzel/dea/stage/Lenna_haar_reconstruction.avi)
|
|
- [CompressionInterp-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1.5-1.avi)
|
|
|
this one shows the image reconstruction from its decomposition on a
|
|
- [CompressionInterp-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-2-1.avi)
|
|
|
*Haar* Basis. Haar basis are the first and simpler example of wavelets.
|
|
- [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)
|
|
|
The following videos are experiments related to the methods described in
|
|
- [CompressionLaplace-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-.5-1.avi)
|
|
|
my report.The two numbers are the corresponding values for alpha and
|
|
- [CompressionLaplace-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-1.avi)
|
|
|
beta.
|
|
- [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)
|
|
|
- [CompressionInterp-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-.5-.5.avi)
|
|
- [CompressionLaplace-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-2.avi)
|
|
|
- [CompressionInterp-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-.5-1.avi)
|
|
- [CompressionLaplace-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1.5-1.avi)
|
|
|
- [CompressionInterp-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-.5.avi)
|
|
- [CompressionLaplace-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-2-1.avi)
|
|
|
- [CompressionInterp-1-1.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-1.5.avi)
|
|
- [CompressionLaplace-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-2-2.avi)
|
|
|
- [CompressionInterp-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-1.avi)
|
|
- [CompressionDisoc-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-2-1.avi)
|
|
|
- [CompressionInterp-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1-2.avi)
|
|
- [CompressionDisoc-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-2-2.avi)
|
|
|
- [CompressionInterp-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-1.5-1.avi)
|
|
- [CompressionDisoc-1.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1.5-1.avi)
|
|
|
- [CompressionInterp-2-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-2-1.avi)
|
|
- [CompressionDisoc-1-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-2.avi)
|
|
|
- [CompressionInterp-2-2.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionInterp-2-2.avi)
|
|
- [CompressionDisoc-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-.5-1.avi)
|
|
|
- [CompressionLaplace-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-.5-.5.avi)
|
|
- [CompressionDisoc-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-1.avi)
|
|
|
- [CompressionLaplace-.5-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-.5-1.avi)
|
|
- [CompressionDisoc-1-1.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-1.5.avi)
|
|
|
- [CompressionLaplace-1-1.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-1.avi)
|
|
- [CompressionDisoc-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-1-.5.avi)
|
|
|
- [CompressionLaplace-1-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionLaplace-1-.5.avi)
|
|
- [CompressionDisoc-.5-.5.avi](http://www.freehackers.org/media/orzel/dea/stage/results/CompressionDisoc-.5-.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)
|
|
## Horror story
|
|
|
- [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)
|
|
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.
|
|
|
- [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)
|
|
## Projet EDP et méthodes mathématiques en finance
|
|
|
- [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)
|
|
 Projet realisé dans le cadre de l'évaluation du cours de Messieurs Pironneau et Berestycki.
|
|
|
- [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)
|
|
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 :
|
|
|
- [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)
|
|
`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.`
|
|
|
- [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)
|
|
|
|
|
|
|
- [Rapport du projet](http://www.freehackers.org/media/orzel/dea/finance/rapport_finance.pdf)
|
|
|
## Horror story
|
|
- [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)
|
|
|
On of the most frightening thing I've ever done with the computer. They
|
|
|
|
|
*requested* me to write something in
|
|
## Projet sur les maillages adaptatifs
|
|
|
[Fortran](http://en.wikipedia.org/wiki/Fortran), so I did an
|
|
|
|
|
[implementation of the
|
|
Projets realises dans le cadre de l'evaluation du cours de Messieurs Hecht et Frey, et réalisé en binôme avec Hanen Amor.
|
|
|
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
|
|
Nous avons choisit de faire deux projets:
|
|
|
programming language anymore.
|
|
|
|
|
|
|
- Le projet 4 : Interpolation de maillages en O(n)
|
|
|
## Projet EDP et méthodes mathématiques en finance
|
|
- 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)
|
|
|
 Projet realisé dans le cadre de
|
|
- [Code source du projet (bz2)](http://www.freehackers.org/media/orzel/dea/maillages/projetmaillage.tar.bz2)
|
|
|
l'évaluation du cours de Messieurs Pironneau et Berestycki.
|
|
- [Code source du projet (tgz)](http://www.freehackers.org/media/orzel/dea/maillages/projetmaillage.tgz)
|
|
|
|
|
|
|
|
J'ai choisi le projet 10, sur les méthodes adaptatives pour la
|
|
### Projet interpolation
|
|
|
résolution de l'équation de Black-Scholes L'énoncé exact est :
|
|
|
|
|
|
|
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é.
|
|
|
` 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 particules
|
|
|
|
|
|
|
|
- [Rapport du
|
|
- [9 particules n'interagissant pas](http://www.freehackers.org/media/orzel/dea/maillages/BouncingMolecule.avi), elles rebondissent juste.
|
|
|
projet](http://www.freehackers.org/media/orzel/dea/finance/rapport_finance.pdf)
|
|
- [validation du tore pour les positions et interactions](http://www.freehackers.org/media/orzel/dea/maillages/ValidationTore.avi) (40 particules)
|
|
|
- [Code source du projet
|
|
- [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)
|
|
|
(bz2)](http://www.freehackers.org/media/orzel/dea/finance/sources-finance.tar.bz2)
|
|
- [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.
|
|
|
- [Code source du projet
|
|
|
|
|
(tar.gz)](http://www.freehackers.org/media/orzel/dea/finance/sources-finance.tar.gz)
|
|
## Projet éléments Finis
|
|
|
|
|
|
|
|
## Projet sur les maillages adaptatifs
|
|
(Cours de M.Perronnet)
|
|
|
|
|
|
|
|
Projets realises dans le cadre de l'evaluation du cours de Messieurs
|
|
L'ensemble du code du projet peut être télécharé ici:
|
|
|
Hecht et Frey, et réalisé en binôme avec Hanen Amor.
|
|
|
|
|
|
|
- [Le projet 2D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem2d.tgz)
|
|
|
Nous avons choisit de faire deux projets:
|
|
- [Le projet 3D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem3d.tgz)
|
|
|
|
|
|
|
|
- Le projet 4 : Interpolation de maillages en O(n)
|
|
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 projet 8 : Schéma numérique et visualisation 3D pour un problème
|
|
|
|
|
d'interaction de particules
|
|
Le deuxième contient l'ensemble du projet dans sa version 3D.
|
|
|
|
|
|
|
|
<!-- -->
|
|
### Premiers essais
|
|
|
|
|
|
|
|
- [Rapport du
|
|
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).
|
|
|
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
|
|
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)
|
|
|
(tgz)](http://www.freehackers.org/media/orzel/dea/maillages/projetmaillage.tgz)
|
|
|
|
|
|
|
|
|
|
### Projet interpolation
|
|
|
|
|
|
|
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.
|
|
|
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
|
|
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)
|
|
|
courant sur le maillage d'entrée. La seconde image montre la même ligne
|
|
|
|
|
sur notre maillage interpolé.
|
|
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:
|
|
|
### Projet particules
|
|
|
|
|
|
|
- [les modifications apportées à Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffFirstTests.html)
|
|
|
- [9 particules n'interagissant
|
|
- [Le fichier complet](http://www.freehackers.org/media/orzel/dea/orzel2d_2trous.cpp)
|
|
|
pas](http://www.freehackers.org/media/orzel/dea/maillages/BouncingMolecule.avi),
|
|
- [Le fichier mesh utilisé](http://www.freehackers.org/media/orzel/dea/orzel_section_2trous.msh)
|
|
|
elles rebondissent juste.
|
|
|
|
|
- [validation du tore pour les positions et
|
|
### Prise en compte du tenseur de conductivité
|
|
|
interactions](http://www.freehackers.org/media/orzel/dea/maillages/ValidationTore.avi)
|
|
|
|
|
(40 particules)
|
|
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.
|
|
|
- [un phenomène intéressant ou des particules s'aggrègent puis
|
|
|
|
|
explosent
|
|
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.
|
|
|
finalement.](http://www.freehackers.org/media/orzel/dea/maillages/agg.avi)
|
|
|
|
|
(40 particules)
|
|
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.
|
|
|
- [une version ou tout
|
|
|
|
|
explose](http://www.freehackers.org/media/orzel/dea/maillages/comportementbizarre.avi),
|
|
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.
|
|
|
pour montrer le type de bug/difficultés que l'on peut rencontrer
|
|
|
|
|
pendant le développement.
|
|
Documents:
|
|
|
|
|
|
|
|
## Projet éléments Finis
|
|
- [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)
|
|
|
(Cours de M.Perronnet)
|
|
|
|
|
|
|
### Conditions de type Neumann
|
|
|
L'ensemble du code du projet peut être télécharé ici:
|
|
|
|
|
|
|
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.
|
|
|
- [Le projet
|
|
|
|
|
2D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem2d.tgz)
|
|
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.
|
|
|
- [Le projet
|
|
|
|
|
3D](http://www.freehackers.org/media/orzel/dea/ThomasCapricelliFem3d.tgz)
|
|
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.
|
|
|
|
|
|
|
|
Le premier contient une copie de fem-v5, légérement modifiée comme
|
|
_Premier Test_
|
|
|
indiqué dans les pages précédentes, les différentes versions du
|
|
|
|
|
programme Laplace, et les fichiers de données correspondants.
|
|
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 :
|
|
|
|
|
|
|
|
Le deuxième contient l'ensemble du projet dans sa version 3D.
|
|
|
|
|
|
|
|
|
|
### Premiers essais
|
|
_Tests_
|
|
|
|
|
|
|
|
J'ai modélisé la coupe transversale d'un tuyau transportant deux fluides
|
|
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:
|
|
|
à 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).
|
|
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
|
|
|
|
|
|
|
|
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
|
|
// test de validité
|
|
|
premier exemple, le gradient est nul sur l'axe y=x, et est linéaire
|
|
R max = 0.;
|
|
|
selon l'axe y=-x, négatif en bas à droite, positif en haut à gauche. (La
|
|
R tetamax =0.;
|
|
|
formule est du type -x+y)
|
|
assert (Ug.N()==Th.nv);
|
|
|
|
|
for ( int i = 0 ; i< Th.nv; i++ ) {
|
|
|
|
|
R diff;
|
|
|
|
|
R2 P;
|
|
|
La température sur l'ensemble de l'élément est majorée par 0 et 10,
|
|
P = (R2)(Th.vertices[i]);
|
|
|
c'est-à-dire les températures extrêmes.
|
|
diff = Ug[i]-TetaExact(P) ;
|
|
|
|
|
diff *= diff; // on prend le carre
|
|
|
J'ai ensuite considéré un gradient nul sur l'axe y=-x, et linéaire selon
|
|
if ( diff>max ) max= diff;
|
|
|
l'axe y=x, négatif en bas à gauche, positif en haut à droite. (La
|
|
if ( Ug[i]>tetamax ) tetamax= Ug[i];
|
|
|
formule est essentiellement x+y)
|
|
}
|
|
|
|
|
printf("Max des carres obtenu : %f\n", max);
|
|
|
Contrairement à l'exemple précédent, la température dépasse les valeurs
|
|
printf("Max des tetas obtenu : %f\n\n", tetamax);
|
|
|
extrêmes, en négatif (bas-gauche) ou en positif (haut-droit).
|
|
```
|
|
|
|
|
|
|
|
|
|
Avec le premier maillage, correspondant à respectivement 16 et 32 segments sur les cercles intérieur et extérieur, j'obtiens :
|
|
|
|
|
|
|
|
Documents:
|
|
`Max des carres obtenu : 0.750861`
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
Avec un maillage plus fin, (50 et 200 segments sur les cercles), j'obtiens :
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffFirstTests.html)
|
|
|
|
|
- [Le fichier
|
|
`Max des carres obtenu : 0.005458`
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/orzel2d_2trous.cpp)
|
|
|
|
|
- [Le fichier mesh
|
|
Ce qui me parait raisonnable pour des températures entre 16 et 100.
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/orzel_section_2trous.msh)
|
|
|
|
|
|
|
Documents:
|
|
|
### Prise en compte du tenseur de conductivité
|
|
|
|
|
|
|
- [les modifications apportées à Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffNeumann.html)
|
|
|
Notre première modification vise à prendre en compte le tenseur de
|
|
- [Le fichier complet](http://www.freehackers.org/media/orzel/dea/LaplaceNeumann.cpp)
|
|
|
conductivité en x et y, c'est-à-dire en fonction du matériau, ou encore,
|
|
- [Le maillage grossier utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere.msh)
|
|
|
en terme Mefisto, du Label.
|
|
- [Le maillage fin utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere_big.msh)
|
|
|
|
|
|
|
|
On s'est contenté ici d'utiliser une matrice de conductivité diagonale,
|
|
### Conditions de type Robin/Fourier
|
|
|
la multiplication par un vecteur revient donc à multiplier chaque
|
|
|
|
|
composante par la composante correspondante de la diagonale.
|
|
L'étape suivante consiste à prendre en compte les conditions de type Robin/Fourier. Pour cela, il faut des données supplémentaires :
|
|
|
|
|
|
|
|
L'équation différentielle vérifiée est alors -div\[A(x)\*grad(t(x))\] =
|
|
- La température extérieure (fluide léchant le bord)
|
|
|
f(x), où t(x) représente la température en x, et A(x) la matrice du
|
|
- Le coefficient d'échange
|
|
|
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
|
|
Dans la formulation variationnelle correspondante, il faut alors modifier à la fois le second membre, et la matrice globale.
|
|
|
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
|
|
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.
|
|
|
noeuds du maillage.
|
|
|
|
|
|
|
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\] <span dir="">\~</span> 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.
|
|
|
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
|
|
_Modifications apportées à fem v5_
|
|
|
deux éléments.
|
|
|
|
|
|
|
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+=)
|
|
|
Documents:
|
|
|
|
|
|
|
```
|
|
|
- [les modifications apportées a
|
|
template<class R> MatriceCreuse<R> & MatriceProfile<R>::add2x2Matrice (int i0, int i1,R a00, R a01, R a11)
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffTenseurConductivite.html)
|
|
{
|
|
|
- [Le fichier
|
|
int i,j;
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceConductivite.cpp)
|
|
if (!D) {
|
|
|
|
|
// matrice vide
|
|
|
### Conditions de type Neumann
|
|
cerr << "Big bug: empty matrice not handled (yet) in add2x2Matrice" << endl;
|
|
|
|
|
exit(1);
|
|
|
Jusqu'à alors, nous ne pouvions utiliser que des conditions de type
|
|
}
|
|
|
Dirichlet. Je vais rajouter, dans un premier temps, la prise en compte
|
|
D[i0] += a00;
|
|
|
de conditions de type Neumann. C'est le plus facile, car il ne faut
|
|
D[i1] += a11;
|
|
|
modifier que le second membre.
|
|
|
|
|
|
|
assert( i0 != i1 );
|
|
|
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
|
|
i=i0; j=i1;
|
|
|
frontière où il y a des conditions de type Neumann. Sinon, bien sûr, je
|
|
if (j<i) L[ pL[i+1] - (i-j) ] += a01;
|
|
|
passe au côté suivant.
|
|
else if (j>i) U[ pU[j+1] - (j-i) ] += a01;
|
|
|
|
|
|
|
|
J'utilise alors une formule d'intégration numérique en 1D pour calculer
|
|
i=i1; j=i0;
|
|
|
l'intégrale sur le segment. Il ne reste plus qu'à modifier le second
|
|
if (j<i) L[ pL[i+1] - (i-j) ] += a01;
|
|
|
membre en conséquence.
|
|
else if (j>i) U[ pU[j+1] - (j-i) ] += a01;
|
|
|
|
|
|
|
|
*Premier Test*
|
|
return *this;
|
|
|
|
|
}
|
|
|
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
|
|
_Tests_
|
|
|
+20 en haut. J'ai imposé une condition de Neumann nulle sur le cercle
|
|
|
|
|
intérieur, je dois donc obtenir des isothermes perpendiculaires à la
|
|
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).
|
|
|
frontière. Ce qui est bien le cas :
|
|
|
|
|
|
|
J'obtiens, avec le maillage fin :
|
|
|
|
|
|
|
|
|
|
`Max des carres obtenu : 0.299266`
|
|
|
*Tests*
|
|
|
|
|
|
|
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.
|
|
|
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
|
|
Documents:
|
|
|
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
|
|
- [les modifications apportées à Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffRobin.html)
|
|
|
téta(x,y)=x\*x+y\*y. Ainsi, l'inverse du laplacien est constant égale à
|
|
- [Le fichier complet](http://www.freehackers.org/media/orzel/dea/LaplaceRobin.cpp)
|
|
|
-4. Je considère alors le problème avec conditions de Dirichlet sur le
|
|
- [Le maillage fin utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere_big.msh)
|
|
|
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
|
|
### Exemple concret en 2D, sans condition de robin
|
|
|
cercle vaut (x,y)/R2, donc la derivée normale vaut 2\*R2). J'obtiens
|
|
|
|
|
ceci:
|
|
Pour traiter un exemple, j'ai choisi de modéliser la température d'une pièce en L dont voici le maillage :
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Ce qui est déjà prometteur, car la température est distribuée
|
|
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.
|
|
|
radialement entre 4\*4 et 10\*10.
|
|
|
|
|
|
|
Bien entendu, l'air de la pièce et les murs n'ont pas la même conductivité.
|
|
|
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.
|
|
Les conditions imposées sont les suivantes :
|
|
|
Voici le code
|
|
|
|
|
|
|
- Température de 0 degré à l'extérieur, sur le mur et la fenêtre de droite (Dirichlet)
|
|
|
` // test de validite`
|
|
- Température de 20 degrés sur les murs intérieurs (Dirichlet)
|
|
|
` R max = 0.;`
|
|
- Le mur et la fenêtre du haut sont baignés par le soleil (Condition de Neumann)
|
|
|
` R tetamax =0.;`
|
|
|
|
|
` assert ( Ug.N() == Th.nv );`
|
|
Voici ce que j'obtiens:
|
|
|
` for ( int i = 0 ; i< Th.nv; i++ )`
|
|
|
|
|
` {`
|
|
|
|
|
` R diff;`
|
|
|
|
|
` R2 P;`
|
|
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.
|
|
|
` `
|
|
|
|
|
` P = (R2)(Th.vertices[i]);`
|
|
_Exemple complet avec condition de robin_
|
|
|
` 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);`
|
|
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)
|
|
|
` printf("Max des tetas obtenu : %f\n\n", tetamax);`
|
|
|
|
|
|
|
Documents:
|
|
|
Avec le premier maillage, correspondant à respectivement 16 et 32
|
|
|
|
|
segments sur les cercles intérieur et extérieur, j'obtiens :
|
|
- [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)
|
|
|
` Max des carres obtenu : 0.750861 `
|
|
- [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)
|
|
|
Avec un maillage plus fin, (50 et 200 segments sur les cercles),
|
|
|
|
|
j'obtiens :
|
|
### Passage à la 3D
|
|
|
|
|
|
|
|
` Max des carres obtenu : 0.005458 `
|
|
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.
|
|
|
|
|
|
|
|
Ce qui me parait raisonnable pour des températures entre 16 et 100.
|
|
_Modifications portant sur le projet_
|
|
|
|
|
|
|
|
Documents:
|
|
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.
|
|
|
|
|
|
|
|
- [les modifications apportées à
|
|
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.
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffNeumann.html)
|
|
|
|
|
- [Le fichier
|
|
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.
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceNeumann.cpp)
|
|
|
|
|
- [Le maillage grossier
|
|
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...
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere.msh)
|
|
|
|
|
- [Le maillage fin
|
|
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.
|
|
|
utilisé](http://www.freehackers.org/media/orzel/dea/testfrontiere_big.msh)
|
|
|
|
|
|
|
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
|
|
|
### Conditions de type Robin/Fourier
|
|
|
|
|
|
|
```plaintext
|
|
|
L'étape suivante consiste à prendre en compte les conditions de type
|
|
R3 Tetraedre::H( int i ) const
|
|
|
Robin/Fourier. Pour cela, il faut des données supplémentaires :
|
|
{
|
|
|
|
|
assert( i>=0 ); assert( i<4 );
|
|
|
- La température extérieure (fluide léchant le bord)
|
|
|
|
|
- Le coefficient d'échange
|
|
// si i est impair, on inverse c et b pour avoir les
|
|
|
|
|
// indices dans l'ordre direct
|
|
|
Dans la formulation variationnelle correspondante, il faut alors
|
|
R3 &D = *(vertices[ i ]),
|
|
|
modifier à la fois le second membre, et la matrice globale.
|
|
&C = *(vertices[ (i+1 + (i%2))%4 ]),
|
|
|
|
|
&B = *(vertices[ (i+2 - (i%2))%4 ]),
|
|
|
Pour le second membre, on fait exactement comme pour les conditions de
|
|
&A = *(vertices[ (i+3)%4 ]);
|
|
|
Neumann, en remplaçant UNeumann(P) par GRobin(P) \* TetaExtRobin(P) avec
|
|
|
|
|
des notations évidentes.
|
|
// check we took the coordinates in the right order
|
|
|
|
|
assert ( ((B-A)^(C-A),(D-A)) >0 );
|
|
|
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
|
|
// On cherche HD, avec H dans le plan ABC et HD perpendiculaire a ce plan
|
|
|
g constant, je peux utiliser la formule intégrale(f) sur \[0,1\] \~
|
|
R3 h = ( B-A )^( C-A ); // normal a la face ABC
|
|
|
1/6(f(0)+4f(1/2)+f(1)), qui est exacte sur les polynômes P2 que je
|
|
R norm = Norme2(h);;
|
|
|
considère ici. Après changement de variable pour se ramener au segment
|
|
h /= norm; // h est unitaire
|
|
|
courant, cela revient à g\*longueur(segment)/6, si i est différent de j,
|
|
h /= (6.* volume/norm); // on divise par la longueur de la hauteur
|
|
|
et à g\*longueur(segment)/3 si i==j.
|
|
return h;
|
|
|
|
|
}
|
|
|
*Modifications apportées à fem v5*
|
|
```
|
|
|
|
|
|
|
|
Pour modifier la matrice globale, je n'ai pas trouvé de fonction
|
|
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..
|
|
|
correspondante dans l'API de fem, j'ai donc rajouté la fonction
|
|
|
|
|
Add2x2Matrice aux fichiers MatriceCreuse_tpl.hpp et MatriceCreuse.hpp,
|
|
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.
|
|
|
dont voici le code (clairement inspiré de la fonction operator+=)
|
|
|
|
|
|
|
_Modifications portant sur le programme principal_
|
|
|
` template`<class R>
|
|
|
|
|
` MatriceCreuse`<R>` & MatriceProfile`<R>`::add2x2Matrice`
|
|
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..
|
|
|
` (int i0, int i1,R a00, R a01, R a11) {`
|
|
|
|
|
` int i,j;`
|
|
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...
|
|
|
` if (!D) { // matrice vide·`
|
|
|
|
|
` cerr << "Big bug empty matrice not handled (yet) in add2x2Matrice" << endl;`
|
|
_Tests_
|
|
|
` exit(1);`
|
|
|
|
|
` }`
|
|
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 <span dir="">#</span>if/#endif dans le code et le rends plus lisible.
|
|
|
` D[i0] += a00;`
|
|
|
|
|
` D[i1] += a11;`
|
|
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.
|
|
|
` `
|
|
|
|
|
` assert( i0 != i1 );`
|
|
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.
|
|
|
` `
|
|
|
|
|
` i=i0; j=i1;`
|
|
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.
|
|
|
` if (j < i) L[ pL[i+1] - (i-j) ] += a01;`
|
|
|
|
|
` else if (j > i) U[ pU[j+1] - (j-i) ] += a01;`
|
|
Enfin le test **donnees_x2y2z2.cpp** est l'exemple classique, et m'a permis de comparer les resultats obtenus avec d'autres eleves.
|
|
|
` `
|
|
|
|
|
` i=i1; j=i0;`
|
|
_Résultats_
|
|
|
` if (j < i) L[ pL[i+1] - (i-j) ] += a01;`
|
|
|
|
|
` else if (j > i) U[ pU[j+1] - (j-i) ] += a01;`
|
|
- **donnees_dirichlet_uniforme.cpp** : j'obtiens exactement la temperature attendue.
|
|
|
` `
|
|
- **donnees_dirichletnul.cpp** : la aussi j'obtiens exactement la temperature attendue.
|
|
|
` return *this;`
|
|
|
|
|
` }`
|
|
`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`
|
|
|
|
|
|
|
|
*Tests*
|
|
- **donnees_quelconque.cpp** : donnees de Dirichlet : La aussi j'obtiens exactement la temperature attendue.
|
|
|
|
|
|
|
|
Pour tester cette partie, j'ai utilisé le même exemple qu'avec les
|
|
`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`
|
|
|
conditions de Neumann. J'impose toujours la température sur la frontière
|
|
|
|
|
intérieure (Dirichlet, donc). Je fixe une température extérieure
|
|
- **donnees_quelconque.cpp** : donnees de Dirichlet, de Neumann et de Robin/Fourier : (c'est le cas le plus difficile)
|
|
|
(TetaR=1000 dans cet exemple). La valeur de g se calcule facilement : g
|
|
|
|
|
= - 2\*R2/(TetaR-R2\*R2).
|
|
`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`
|
|
|
|
|
|
|
|
J'obtiens, avec le maillage fin :
|
|
### Exemple 3D
|
|
|
|
|
|
|
|
` Max des carres obtenu : 0.299266`
|
|
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).
|
|
|
|
|
|
|
|
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:
|
|
|
|
|
|
|
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".
|
|
|
- [les modifications apportées à
|
|
|
|
|
Laplace.cpp](http://www.freehackers.org/media/orzel/dea/DiffRobin.html)
|
|
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.
|
|
|
- [Le fichier
|
|
|
|
|
complet](http://www.freehackers.org/media/orzel/dea/LaplaceRobin.cpp)
|
|
**UPDATE :** M.Perronnet m'a expliqué le problème : la surface inférieure de l'hexaèdre de la résistance est confondu avec le disque inférieur du cylindre, et cela pose problème à Mefisto. A la rigueur, si le maillage correspondait, cela serait peut-être passé... |
|
|
- [Le maillage fin
|
|
\ No newline at end of file |
|
|
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é... |
|
|
