sur une approche à objets généralisée pour la mécanique non linéaire

UNIVERSITÉ DE PROVENCE (AIX-MARSEILLE I)

Thèse de Doctorat

Discipline : Mécanique

SUR UNE APPROCHE À OBJETS GÉNÉRALISÉE

POUR LA MÉCANIQUE NON LINÉAIRE

présentée et soutenue publiquement

par

ROY SAAD

le 05 Décembre 2011

dirigée par :

Dominique Eyheramendy

JURY

M. D. Caromel Professeur, INRIA, Université de Nice Examinateur

M. J-Y. Cognard Professeur, ENSTA de Bretagne Rapporteur

M. D. Eyheramendy Professeur, École Centrale de Marseille Directeur

M. F. Feyel Maître de recherches, ONERA, École des Mines de Paris Rapporteur

M. F. Lebon Professeur, Université d‟Aix Marseille Président

M. J. Liandrat Professeur, École Centrale de Marseille Examinateur

M. C. Rey Professeur, École Normale Supérieure de Cachan Examinateur

A ma famille…

A mes amis…

Remerciements

Remerciements

Je tiens tout d’abord à remercier Dieu pour la force et la patience qu’il m’a donné pour

mener à terme ce travail.

Cette thèse a été réalisée au Laboratoire de Mécanique et d’Acoustique à l’École Centrale de

Marseille, sous la direction de Dominique Eyheramendy que je souhaite remercier pour l'aide

compétente qu'il m'a apportée, ses conseils bénéfiques, sa grande disponibilité et son soutien

qui m’a été bien précieux. J'aimerais également lui dire à quel point j’ai apprécié ses qualités

professionnelles et humaines tout au long de ce travail.

Je remercie sincèrement les membres de mon jury : Messieurs Jean-Yves Cognard et Frédéric

Feyel pour avoir accepté de rapporter mon mémoire et Messieurs Christian Rey, Frédéric

Lebon, Jacques Liandrat et Denis Caromel pour avoir accepté de participer à ce jury.

Je souhaite également adresser mes remerciements à l’ensemble des membres de l’équipe

Modèles Numériques du Laboratoire de Mécanique et d’Acoustique : Adnane, Bruno,

Thierry, Stéphane, Jean et Stéphane pour leurs conseils tout au long de ces années. Je ne

saurais oublier l’ensemble des thésards que j’ai rencontrés. Ceux avec lesquels j’ai eu la

chance de partager mon bureau : François, Elia et surtout Thienan pour les différents

échanges qui m’ont été très utiles.

Je n’oublierai pas le groupe de midi (et de soir aussi) : Georges, Kifah, Mansour, Marjorie,

Moncef, Audrey et Mylène. Merci pour votre amitié et pour tous les moments passés

enembles.

Enfin, mes plus vifs remerciements s’adressent à ma mère Rosalie, mon père Miled, et mes

deux frères Abdo et Ralph pour leur soutien et leur encouragement au cours de ces années.

Résumé

Résumé

Les problèmes qui se posent aujourd'hui en mécanique numérique et domaines connexes sont

complexes, et impliquent de plus en plus souvent plusieurs physiques à différentes échelles de

temps et d‟espace. Leur traitement numérique est en général long et difficile, d‟où l‟intérêt

d‟avoir accès à des méthodes et outils facilitant l‟intégration de nouveaux modèles physiques

dans des outils de simulation. Ce travail se pose dans la problématique du développement de

codes de calcul numérique. L‟approche proposée couvre la démarche de développement du

modèle numérique depuis la formulation variationnelle jusqu‟à l‟outil de simulation.

L‟approche est appliquée à la méthode des éléments finis. Nous avons développé des concepts

génériques afin d‟automatiser la méthode des éléments finis. Nous nous sommes appuyés sur

l'analyse tensorielle dans le contexte de la méthode des éléments finis. Le formalisme

mathématique est basé sur l‟algèbre tensorielle appliquée à la description de la discrétisation

des formes variationnelles. Ce caractère générique est conservé grâce à l'approche logicielle

choisie pour l‟implantation; orientée objet en Java. Nous proposons donc un cadre orienté

objet, basé sur des concepts symboliques, capables de gérer de manière symbolique les

développements assistés des contributions élémentaires pour la méthode éléments finis. Ces

contributions sont ensuite automatiquement programmées dans un code de calcul. L'intérêt de

cette approche est la généricité de la description qui peut être étendue naturellement à tout

autre modèle de discrétisation (spatiale ou temporelle). Dans ce travail, les concepts sont

validés dans le cadre de problèmes linéaires simples (élasticité, chaleur,...), dans le cadre du

traitement de formulations variationnelles mixtes (thermomécanique, Navier-Stokes,…) et

dans un cadre Lagrangien (élasticité en grandes transformations, hyperélasticité,…).

Abstract

Abstract

The problems occurring today in computational mechanics and related domains are complex,

and may involve several physics at different time and space scales. The numerical treatment

of complex problems is in general tough and time consuming. In this context, the interest to

develop methods and tools to accelerate the integration of new formulations into simulation

tools is obvious. This work arises on the issue of the development of computational tool. The

proposed approach covers the development process of numerical models from the variational

statement to the simulation tool. The approach is applied to the finite element method. We

have developed generic concepts to automate the development of the finite element

method. To achieve this goal, we relied on tensor analysis applied in the context of the finite

element method. The mathematical formalism is based on the tensor algebra to describe the

discretization of a variational formulation. The generic character of the approach is preserved

through the object-oriented approach in Java. We propose a framework based on object-

oriented concepts capable of handling symbolic developments of elemental contributions for

finite element codes. The advantage of this approach is the generic description that can be

extended naturally to any discretization model in space or time. This concept is fully validated

for simple linear problems (elasticity, heat convection, ...), for the treatment of mixed

variational formulations (thermo-mechanical, Navier-Stokes for incompressible flows...) and

Lagrangian frameworks (elasticity in larges transformations, hyperelasticity, ...).

Tables des matières

11

Table des matières

CHAPITRE 1 INTRODUCTION .................................................................. 17

CHAPITRE 2 LANGAGES A OBJETS, CODES ELEMENTS FINIS ET

APPROCHES SYMBOLIQUES ..................................................................... 21

2.1 L’EVOLUTION DES LANGAGES DE PROGRAMMATION POUR L’IMPLEMENTATION DE LA

METHODE DES ELEMENTS FINIS ....................................................................................... 21

2.2 PROGRAMMATION ORIENTEE OBJET ET CODE ELEMENTS FINIS .................................... 23

2.2.1 Les précurseurs en matière de code orientés objet en mécanique .................... 23

2.2.2 La maturité des langages orientés objet appliqués à la mécanique .................. 24

2.2.3 Vers des approches structurées plus avancées .................................................. 24

2.3 L’UTILISATION DU CALCUL SYMBOLIQUE DANS LES APPLICATIONS ELEMENTS FINIS .... 26

2.3.1 Approches semi-analytiques numériques ......................................................... 26

2.3.2 Amélioration de la performance des codes éléments finis ............................... 27

2.3.3 Efficacité et flexibilité dans le développement de codes éléments finis........... 27

2.4 UNE APPROCHE ORIENTEE OBJET INTEGREE SYMBOLIQUE/NUMERIQUE POUR

L’ANALYSE ELEMENTS FINIS ............................................................................................. 31

CHAPITRE 3 UNE APPROCHE TENSORIELLE GENERIQUE POUR

LA DERIVATION DES MODELES ELEMENTS FINIS ........................... 35

3.1 DESCRIPTION MATHEMATIQUE DE L’APPROCHE DE DISCRETISATION ELEMENTS

FINIS…… ....................................................................................................................... 36

3.1.1 Formulation variationnelle ................................................................................ 36

3.1.2 Discrétisation éléments finis ............................................................................. 36


12

3.1.3 Représentation tensorielle ................................................................................. 37

3.2 ILLUSTRATION DU FORMALISME SUR DES FORMULATIONS SIMPLES ............................. 38

3.2.1 Application à l‟élasticité linéaire ...................................................................... 38

3.2.2 Application à un problème non linéaire simple ................................................ 41

3.3 ANALYSE PRELIMINAIRE SUR L’IMPLANTATION LOGICIELLE DE LA DISCRETISATION

DES FORMES VARIATIONNELLES ...................................................................................... 47

CHAPITRE 4 UN ENVIRONNEMENT ORIENTE OBJET POUR LA

DERIVATION SYMBOLIQUE DE FORMULATIONS ELEMENTS

FINIS…………… .............................................................................................. 49

4.1 UN ENVIRONNEMENT A OBJETS POUR LE CALCUL FORMEL ............................................ 49

4.2 ANALYSE A OBJET POUR LA REPRESENTATION SYMBOLIQUE D’UNE FORME

VARIATIONNELLE ............................................................................................................. 51

4.3 DISCRETISATION DE LA FORME VARIATIONNELLE ET FORMES ELEMENTAIRES............ 65

4.3.1 Discrétisation symbolique ................................................................................ 65

4.3.2 Génération des contributions élémentaires ....................................................... 68

CHAPITRE 5 GENERATION AUTOMATIQUE DE CODE

ELEMENTS FINIS ........................................................................................... 73

5.1 GENERATION AUTOMATIQUE DE CODE ET CODE ELEMENTS FINIS ................................. 73

5.2 UN CODE ELEMENTS FINIS ORIENTE OBJET EN JAVA : FEMJAVA ................................. 74

5.2.1 Organisation des packages ................................................................................ 75

5.2.2 Développement de nouveaux modèles éléments finis dans FEMJava ............. 79

5.2.2.1 Principe d‟intégration d‟un nouvel élément dans FEMJava .............. 79

5.2.2.2 Application à un écoulement d‟un fluide incompressible .................. 80

5.3 OUTIL DE GENERATION AUTOMATIQUE DE CODE ........................................................... 83

5.3.1 Développement d‟une formulation variationnelle sous forme symbolique ...... 83

5.3.1.1 Définition des éléments de base du problème .................................... 85

5.3.1.2 Construction des termes de base du problème ................................... 88


13

5.3.1.3 Discrétisation des champs solutions et tests ...................................... 88

5.3.1.4 Les opérations tensorielles pour construire les matrices

élémentaires ........................................................................................ 89

5.3.1.5 Définition de la formulation du problème .......................................... 90

5.3.2 Génération automatique de code à partir des formes symboliques .................. 90

5.3.2.1 Intégration de l‟environnement symbolique dans l‟environnement

numérique .......................................................................................... 93

5.3.2.2 Particularisation les champs inconnus au problème .......................... 93

5.3.2.3 Entêtes de méthodes pour le calcul des grandeurs élémentaires ........ 94

5.3.2.4 Calcul des contributions élémentaires aux points d‟intégration ........ 95

5.3.3 Application numérique ..................................................................................... 97

CHAPITRE 6 APPROCHE SYMBOLIQUE POUR FORMES

VARIATIONNELLES EN GRANDES DEFORMATIONS ........................ 99

6.1 APPROCHES SYMBOLIQUES ET GRANDES TRANSFORMATIONS ...................................... 99

6.2 CINEMATIQUE DES GRANDES TRANSFORMATIONS ....................................................... 100

6.3 FORMULATIONS FORTE ET FAIBLE D’UN PROBLEME EN GRANDES

TRANSFORMATIONS ....................................................................................................... 101

6.3.1 Formulation forte ............................................................................................ 101

6.3.2 Formulation faible .......................................................................................... 101

6.3.3 Formulation faible linéarisée .......................................................................... 102

6.4 GENERATION DES CONTRIBUTIONS ELEMENTAIRES ..................................................... 104

6.5 DESCRIPTION SYMBOLIQUE DU PROBLEME ................................................................... 107

6.5.1 Définition des éléments de base de la formulation ......................................... 108

6.5.2 Définitions des éléments discrétisés et contributions élémentaires ................ 109

6.5.3 Génération de code ......................................................................................... 112

6.5.4 Application numérique ................................................................................... 114

CHAPITRE 7 APPLICATIONS ................................................................. 117

7.1 HYPERELASTICITE ......................................................................................................... 117


14

7.1.1 Formulation Forte ........................................................................................... 117

7.1.2 Formulation faible .......................................................................................... 118

7.1.3 Linéarisation de la forme faible ...................................................................... 119

7.1.4 Description symbolique du problème ............................................................. 120

7.1.4.1 Définition des éléments de base de la formulation .......................... 120

7.1.4.2 Définitions des éléments discrétisés et contributions élémentaires . 122

7.1.5 Application numérique ................................................................................... 123

7.2 THERMOELASTICITE ...................................................................................................... 124

7.2.1 Loi de comportement mécanique et loi de Fourier ......................................... 124

7.2.2 Equation de la chaleur .................................................................................... 127

7.2.3 Forme forte et forme faible ............................................................................. 129

7.2.4 Description symbolique du problème ............................................................. 129

7.2.4.1 Définition des éléments de base de la formulation .......................... 129

7.2.4.2 Définitions des éléments discrétisés et contributions élémentaires

associées .......................................................................................... 131

7.2.5 Applications numériques ................................................................................ 132

CHAPITRE 8 CONCLUSION .................................................................... 137

8.1 SUR UNE APPROCHE A OBJETS GENERALISEE POUR LA MECANIQUE NON LINEAIRE .... 137

8.1.1 Cadre général de l‟approche ........................................................................... 137

8.1.2 Une approche tensorielle générique pour la dérivation de modèles éléments

finis ................................................................................................................. 138

8.1.3 Un environnement à objets de calcul formel pour la dérivation de formulations

éléments finis .................................................................................................. 139

8.1.4 Une approche de génération automatique de code éléments finis .................. 139

8.1.5 Sur la validation de l‟approche proposée ........................................................ 140

8.2 POUR ALLER PLUS LOIN DANS CE TRAVAIL ................................................................... 141

8.3 LES OUTILS DE CALCUL DE DEMAIN .............................................................................. 143


15

ANNEXE 1 JAVA, UN LANGAGE A OBJETS POUR LE CALCUL

SCIENTIFIQUE .............................................................................................. 145

ANNEXE 2 DERIVATION DES EQUATIONS DE NAVIER-STOKES EN

MILIEU INCOMPRESSIBLE....................................................................... 159


16

Chapitre 1 Introduction

17


En ingénierie moderne, la simulation numérique est devenue aujourd'hui un outil

incontournable. Les problèmes auxquels les ingénieurs et les chercheurs s'attaquent

aujourd'hui deviennent de plus en plus complexes. Une des tendances fortes est de développer

des modèles et outils de simulation faisant intervenir différentes physiques à la mécanique. Le

formalisme permettant de décrire les grandes transformations ajoutent un niveau de

complexité supplémentaire. Ces outils de simulation sont intégrés dans des systèmes

informatiques pouvant eux-mêmes être complexes (systèmes multiprocesseurs, grilles de

calcul,...). Un ingénieur ou chercheur développant un nouveau modèle en mécanique, incluant

par exemple des physiques complexes couplées, est confronté au choix de l'outil numérique

afin d‟y intégrer son modèle. Suivant la taille du modèle (nombre de degrés de libertés,

nombre de champs couplés, contraintes, déplacement, température, pression, densité, taux

d‟humidité,…) et les couplages entre les différentes physiques, il est nécessaire de développer

ou d'adapter les algorithmes et méthodes numériques de résolution existantes.

Dans ce contexte, le développement de code éléments finis s‟avère être une activité difficile

dés lors que l‟on traite des formulations complexes. Le point de départ de la démarche est un

problème physique. Ce problème est généralement modélisé par un ensemble d‟équations :

équations aux dérivées partielles, équation aux dérivées ordinaires,… Une procédure éléments

finis est alors appliquée au modèle mathématique. La démarche traditionnelle consiste alors à


18

élaborer un outil informatisé, en général complexe et parfois relativement éloigné des

descriptions mathématiques. En général, il est laissé à la charge de l‟utilisateur final le soin de

modéliser son problème particulier: géométrie, contraintes cinématiques, charges, données du

problème… Dans ce contexte, les approches structurées de programmation, et en particulier à

objets, ont permis de faire évoluer très notablement les outils de simulation en mécanique.

D'un point de vue technique, vu la complexité des modèles étudiés, on cherche souvent à

s'appuyer sur des outils existants (Abaqus, Nastran,...) qui bien souvent freinent le

développement d'outils numériques, n‟étant pas forcement adaptés aux nouvelles

formulations. Une démarche qui commence à être aujourd‟hui de plus en plus utilisée et qui

permet d‟aller plus loin en général dans le cadre de problèmes couplés, est de s‟appuyer sur

des outils de plus haut niveau d‟abstraction tel que Comsol. Dans ce cas, on se trouve bien

souvent limité par l‟algorithmique globale du solveur. Il est clair au regard de toutes ces

approches possibles que le besoin d‟outils de haut niveau d‟abstraction doit aider au

développement d‟outils de simulation.

Ce travail se pose dans la problématique de développement de codes de simulation.

L‟approche proposée couvre la démarche de développement du modèle numérique depuis la

formulation variationnelle jusqu‟à l‟outil de simulation. L‟approche est appliquée à la

méthode des éléments finis. Nous avons développé des concepts génériques pour automatiser

la méthode des éléments finis. Le formalisme mathématique se base sur l‟algèbre tensorielle

pour décrire la discrétisation d‟une formulation variationnelle. Ce caractère générique est

conservé grâce à l'approche logicielle choisie, orientée objet en Java. Nous proposons donc un

cadre orienté objet intégrant des concepts symboliques capables de gérer les développements

assistés de codes éléments finis. L'approche se décompose en deux étapes. La première étape

permet la construction de la forme variationnelle, puis la discrétisation des champs adéquats

et enfin, la génération sous forme symbolique des contributions élémentaires basé sur

l‟algèbre tensorielle. La deuxième étape est la pré-compilation automatique de ces formes

élémentaires afin de les intégrer dans un code de calcul classique.

Dans le Chapitre 2, nous présentons d‟un point de vue général la problématique du

développement de code de calcul en mécanique ou disciplines connexes. Une étude

bibliographique décrit les approches logicielles de structuration des codes de calcul, en

particulier orientées objets. Les approches de plus haut niveau d‟abstraction, approches

symboliques, pour la dérivation des modèles éléments finis sont aussi décrites. Enfin, nous


19

proposons une approche orientée objet intégrée, symbolique/numérique, pour l‟analyse

éléments finis.

Dans le Chapitre 3, nous décrivons les concepts d‟approche générique pour la dérivation de

formulations éléments finis. Il s'agit de poser les bases d'un environnement logiciel permettant

le développement automatisé d'outils de simulation basés sur la dérivation de formes

variationnelles dans un contexte éléments finis. Nous intégrons dans l‟approche tous les

aspects nécessaires pour le traitement de problèmes couplés en mécanique. Nous décrivons les

différentes étapes de l‟approche. Nous présentons le formalisme mathématique sur lequel

nous développerons les outils symboliques incluant ce cadre formel. Puis, nous illustrons le

formalisme mathématique sur deux exemples simples d‟équations aux dérivées partielles :

l‟élasticité linéaire et une équation non linéaire.

Dans le Chapitre 4, nous développons un environnement capable de :

- représenter des formes variationnelles,

- représenter des contributions élémentaires et de les gérer dans le contexte de la

formulation,

- réaliser des calculs formels pour obtenir ces contributions élémentaires

- et enfin générer un code de calcul élément finis.

Cette approche s‟appuie sur le développement d‟un moteur de calcul formel pour l‟algèbre

tensorielle, permettant d‟obtenir sous forme symbolique toutes les contributions élémentaires

d‟un problème dérivé dans un cadre variationnel. Un des points clés de l‟approche que nous

proposons est de n‟utiliser qu‟un seul et même environnement de programmation, orienté

objets en Java. Les objets fondamentaux nécessaires à la représentation symbolique d‟une

forme variationnelle sont présentés en s‟appuyant sur l‟élasticité linéaire ainsi que la

hiérarchie des classes établie de cette étude. Les objets nécessaires à l‟automatisation de la

discrétisation sont aussi présentés.

Dans le Chapitre 5, les concepts essentiels de la génération automatique de code éléments

finis sont donnés. Les contributions élémentaires sont intégrées dans le code éléments finis

FEMJava. Nous décrivons très brièvement le code FEMJava en présentant les principaux

packages du code (développé en Java). Puis nous présentons l‟implémentation de l‟outil de

génération automatique de code. La description s‟appuie sur l‟élasticité linéaire.


20

Dans le Chapitre 6, nous rappelons brièvement le formalisme Lagrangien permettant de

traiter les grandes transformations en mécanique. Nous mettons en évidence les concepts de

base afin de manipuler les formes variationnelles en grandes transformations à partir d‟un

exemple d‟une formulation classique issue de la littérature.

Afin de montrer la généricité de l‟approche proposée dans le Chapitre 7, nous appliquons

l‟approche à deux problèmes de la mécanique : l‟hyperélasticité et la thermoélasticité. Après

avoir présenté les modèles, nous dérivons le modèle numérique dans l‟environnement

symbolique et générons le code éléments finis. Des tests numériques simples valident

l‟approche.

Chapitre 2 Langages à objets, codes éléments finis et approches symboliques

21

Chapitre 2 Langages à objets, codes

éléments finis et approches symboliques

Le développement de code éléments finis s‟avère être une activité difficile dès que l‟on traite

des formulations complexes qui ne s‟intègrent pas naturellement dans des outils prévus à cet

effet. Le point de départ de la démarche est un problème physique. Ce problème est

généralement modélisé par un ensemble d‟équations : équations aux dérivées partielles,

équation aux dérivées ordinaires,… Une procédure éléments finis est alors appliquée au

modèle mathématique. La démarche traditionnelle consiste alors à élaborer un outil

informatisé, en général complexe et parfois relativement éloigné des descriptions

mathématiques. En général, il est laissé à la charge de l‟utilisateur final le soin de modéliser

son problème particulier: géométrie, contraintes cinématiques, charges, données du

problème… Dans ce contexte, les approches structurées de programmation et en particulier à

objets ont permis de faire évoluer très notablement les outils de simulation en mécanique.

2.1 L’évolution des langages de programmation pour

l’implémentation de la méthode des éléments finis

Les langages de programmation généraux ont connu une évolution importante depuis les

années 50. Le premier langage que l‟on peut qualifier de « haut niveau d‟abstraction » a été

Fortran. Ses instructions restent proches dans la syntaxe des formules mathématiques. A sa


22

création, le langage était bien adapté aux besoins des scientifiques. Des langages permettant

une programmation plus structurée furent introduits dans les années 70, parmi lesquels on

peut citer Pascal (1970). C‟est un langage structuré conçu de façon à éviter les erreurs de

programmation en encourageant la modularité. Ce langage fut largement utilisé par les

pédagogues voulant donner aux étudiants une formation initiale à la programmation. Les

langages orientés objet tels que Simula (1969), Smalltalk (1980), C++ (1983) permettent

d'écrire des logiciels fondés sur des objets réutilisables, entités rassemblant données et

traitements, et communiquant entre eux par envoi de messages. En temps, l‟informatique à

vocation scientifique a suivi une évolution parallèle, le tournant se situant avec l‟avènement

des langages à objets dans les années 90.

Le langage traditionnellement utilisé dans la communauté scientifique a longtemps été

FORTRAN. Ce langage permet à l‟origine une programmation modulaire très efficace grâce

au développement des compilateurs associés aux systèmes d‟exploitation, mais ne permettant

aucune utilisation de structures de données de haut niveau pour les versions antérieures à

FORTRAN 95. Dans le domaine de la mécanique, la complexité des programmes augmente

rapidement avec le type de problèmes traités, ainsi qu‟avec la complexité des problèmes

abordés. La maintenance et surtout l‟extension des codes deviennent dans ce contexte très

difficile, un recodage de ces systèmes n'enlevant pas ce manque de flexibilité à moyen terme.

Dès les années 90, le concept orientée objet a permis de gérer cette problématique dans sa

globalité. Le concept clé de la programmation orientée objet est celui d‟objet. Il s‟agit d‟une

entité qui contient à la fois des données (attributs) et des actions (méthodes). Les méthodes

définissent les actions à réaliser au sein d‟un objet. Cela correspond au code exécuté lorsque

l‟objet reçoit un message envoyé par un autre objet. La communication entre objets se fait par

envoi des messages. Les messages que peut recevoir un objet définissent son interface. Cela

définit les seules opérations que l‟on peut demander à un objet de réaliser. Cela empêche

l‟utilisateur d‟accéder directement aux détails de programmation de l‟objet (attributs et

méthodes), qui peuvent évoluer au cours du temps lors de l‟exécution du programme. Cela

s‟appelle l‟encapsulation des données. Les objets ne sont alors que des instances de classes.

Les classes sont organisées en une hiérarchie qui permet l‟héritage. On peut ainsi par héritage

créer une nouvelle classe à partir d'une classe existante. Les attributs et les méthodes ajoutés à

la classe dérivée viennent alors s'ajouter à ceux hérités de la superclasse.


23

L‟utilisation d‟outils d‟informatiques de haut niveau d‟abstraction est indispensable au

développement d‟outils de simulation en mécanique comme montré dans [Fritzon 92]. Ces

outils peuvent être regroupés en trois classes. La première correspond à la plus récente

génération d‟outils de programmation de haut niveau que l‟on peut décomposer en deux

groupes : les langages de programmation procéduraux classiques (FORTRAN 77, 90,

PASCAL,…) et les langages orientés objets (SMALLTALK, C++, Java,…). La deuxième

classe regroupe tous les logiciels mathématiques de manipulations algébriques tels que Maple,

Mathematica, Matlab,…. Enfin, l‟approche hybride pour la résolution des problèmes en

mécanique, ce qui signifie une approche basée sur les outils numériques symboliques mixtes.

2.2 Programmation orientée objet et code éléments finis

2.2.1 Les précurseurs en matière de code orientés objet en mécanique

L‟utilisation de l‟approche orientée objet dans la programmation de codes de simulation par

éléments finis est apparue à la fin des années 80. Elle a été évoquée par [Miller 88] et

appliquée dans [Rehak, 89]. Dans ces deux papiers, les concepts de base pour le calcul de

structure et implémentation dans des systèmes LISP sont donnés. Dans [Fenves 90], les

propriétés de modularité et de réutilisation de code sont mises en avant. Dans [Forde 90],

l‟efficacité dans la maintenance et dans l‟implantation est décrite comme point clé. Les

programmes développés restent cependant limités à des formulations de l‟élasticité linéaire.

Des approches complètes ont été développées pour l‟analyse statique et dynamique des

éléments finis [Zimmerman 92], [Baugh 92], [Scholz 92], [Devloo 94] et [Dubois-Perlin 93].

Dans [Zimmerman 92] et [Dubois-Perlin, 93] la convivialité du Smalltalk pour la phase du

développement mais aussi son manque d‟efficacité par rapport au C++ sont montrés. Par la

suite, des programmes plus complexes ont été développés par [Foerch 95] [Besson 97] (code

Zebulon qui décrit l‟objet matériau de manière avancée et géométrique ainsi que le contact

entre solides avec frottement). On peut tout de même noter que dans le même temps, des

démarches pouvant s‟apparenter à une démarche orientée objet ont été adoptées tout en

gardant FORTRAN comme langage de programmation (par exemple pour les programmes

SIC (Système Interactif de Conception) [Golay 02] et CAST3M, anciennement

CASTEM2000 [Verpeaux 88]). Des modules ont été écrits pour reproduire des fonctionnalités

avancées à la base des langages orientés objets (allocation dynamique de la mémoire, notion

de structure de données,..).


24

2.2.2 La maturité des langages orientés objet appliqués à la mécanique

Depuis cette origine, la programmation orientée objet a été largement appliquée à tous les

domaines du calcul mécanique (voir par exemple [Mackerle 04] pour une vision assez

exhaustive du domaine dans les années 2000). On pourrait à ce stade citer un grand nombre de

travaux. A titre d‟exemple, nous donnerons : problèmes elasto-plastiques [Menetrey 93],

outils numériques comme une bibliothèque d'algèbre linéaire [Zeglinski 94], une version

orientée objet (LAPACK++) de la bibliothèque LAPACK [Dongarra 94], création de codes

interactifs avec l'interface graphique d'utilisateur (GUI) [Mackie 00], l‟intégration de

l'intelligence artificielle dans les systèmes d'éléments finis, en utilisant une architecture qui

partage les données [Bomme et Zimmerman 1996], problèmes de fractures et

d‟endommagement [Hang 02], calcul parallèle en mécanique des solides et fluides [Adeli 00],

mécanique des solides déformables en grandes transformations (Herezh++) [Rio 08] et pour

des problèmes multiphysiques [Eyheramendy 06] et [Wang 10], problèmes de contacts [Fang

00] et [Ma 02]. Cette liste est loin d‟être exhaustive, mais montre que le paradigme objet est

aujourd‟hui largement répandu dans la communauté scientifique du calcul. Dans la majeure

partie des développements, le langage C++ est utilisé surtout pour son efficacité. Mais, on

peut trouver aussi d‟autres approches.

2.2.3 Vers des approches structurées plus avancées

Les langages de programmation de haut niveau d‟abstraction (C, Fortran,…) ont globalement

bien répondu aux exigences de la simulation et du calcul numérique. Par contre, lorsqu‟il

s‟agit de gérer des applications de plus grandes tailles, dans un contexte multi-développeurs

couvrant des domaines plus variés dans le cadre de la mécanique ou autres, des difficultés de

gestion de code apparaissent : maintenance, extension,… Parmi les caractéristiques qui

peuvent permettre d‟améliorer la gestion des applications informatiques, on peut trouver : des

structures de langage supplémentaires bien identifiées telles que les interfaces en Java, des

systèmes de gestion automatique de mémoire (Smalltalk, Java,…), l‟accès à des bibliothèques

avancées (Interfaces Graphiques d'Utilisateur, communication avec d'autres ordinateurs sur un

réseau, internet,…). On peut citer parmi les dernières générations de langages qui sont utilisés

en informatique scientifique, Java et C# ou des langages plus anciens tels que Smalltalk.

Le langage de programmation, Java créé par Sun Microsystems, possède des caractéristiques

qui le rendent attrayant pour écrire des applications de calcul numérique en mécanique. Il est

un langage de programmation simple (plus simple que C ++). Il possède une collection riche

de bibliothèques. Un des principaux avantages de ce langage est sa portabilité. Une


25

application développée sur un système peut être directement exécutée sur toutes les

plateformes, à la restriction de version près. La portabilité de code est obtenue grâce aux

Machines Virtuelles Java (JVM). Les machines virtuelles Java sont développées pour tous les

systèmes informatiques importants. On les trouve également implantées dans les navigateurs

Web ce qui permet d‟exécuter du code Java sur toute plateforme ayant un navigateur ainsi

configuré. Les ''applets'' peuvent être téléchargés par internet et exécutés dans le navigateur de

Web.

Le langage Java est un langage orienté objet. Les notions d‟objet, de classe et héritage y sont

donc implémentées. On y trouve un certain nombre de concepts supplémentaires permettant

d‟améliorer considérablement la structuration des codes dans le cadre de codes de grandes

tailles développés par des équipes de taille importante. Enfin, Java possède suivant les types

de plateforme, un grand nombre de classes prédéfinies au sein de paquetages thématiques. Il

est par exemple facile de créer des Interfaces Graphiques d'Utilisateur (GUI) et de créer des

applications communicantes sur réseau. Le langage Java est présenté de manière succincte

dans l‟Annexe 1 (pour plus de détails voir par exemple [Flanagan 04] et [Flanagan 06]).

Le langage C# développé par Microsoft est très similaire à Java. Il est intégré à la plate-forme

«Microsoft.Net» qui permet au développeur d‟interagir avec d‟autres langages de

programmation (intégration de code développé dans plusieurs langages) ; le CLR (Common

Langage Runtime) facilite les communications y compris dans un contexte distribué. Il reste

cependant bien moins utilisé que le langage Java dans le contexte du calcul scientifique (voir

par exemple [Mackie 01]).

Récemment, des auteurs ont proposé différentes implémentations de la méthode des éléments

finis dans différents domaines de la mécanique en Java et C#. Le projet Carta Blanca

[Vanderheyden et Padial-Collins 03 et 04], outils de simulation de problèmes physiques non

linéaires couplés, en mécanique des fluides [Riley, 03], en mécanique de la rupture

[Nikishkov 99], en calcul parallèle de dynamique de fluide et en mécanique des solides

[Eyheramendy 04, 06, 07], [Baduel 04] en électromagnétisme, [Hauser 00] pour le calcul

parallèle en mécanique [Mackie et Heng 06]. Dans ce dernier travail, l‟implémentation est

réalisée en C# (voir aussi [Mackie 07]); cela représente la seule référence d‟un travail réalisé

en C#. Notons que de nombreux auteurs ont préconisé l'utilisation de Java et ont comparé ses

performances à d‟autres langages de programmation (langage C): [Nikishkov 03, 06 et 10],

[Bull 01], [Eyheramendy 03, 08, 11].


26

2.3 L’utilisation du calcul symbolique dans les applications

éléments finis

L‟utilisation de logiciels de manipulations algébriques a toujours été un centre d‟intérêt pour

le développement d‟élément finis. Les premiers travaux mentionnant leur utilisation datent du

début de l‟utilisation intensive de la méthode des éléments finis dans les années 70. Parmi les

premiers travaux signalés, on trouve [Luft 71]. Dans cet article, une méthodologie pour

générer automatiquement des formes matricielles pour la méthode éléments finis est décrite

en introduisant simplement les caractéristiques du nouvel élément.

Depuis, de nombreux auteurs ont utilisé de telles approches dans le but d‟aider au

développement de procédures de solutions éléments finis. On peut grouper les approches

symboliques dans trois catégories [Eyheramendy 98]. Dans la première, sont regroupés les

travaux dans lesquels les auteurs ont développé des procédures de résolution éléments finis

dans des environnements de calcul algébrique. La deuxième correspond aux approches dont le

but principal est d‟améliorer l‟efficacité numérique des codes éléments finis classiques. Enfin,

de nombreux auteurs ont cherché à accélérer le développement de code de calcul numérique à

l‟aide d‟outils de manipulations symboliques.

2.3.1 Approches semi-analytiques numériques

Dans cette approche, une approche éléments finis classique est utilisée au sein d‟un

environnement de calcul symbolique classique comme Maple, Mathematica, Matlab, etc.

Ainsi, certaines variables peuvent être conservées sous forme symbolique, ce qui permet

d‟évaluer leur influence sur les résultats numériques. Dans [Choi 92], une application pour

l‟élasticité 2D est développée dans l‟environnement Mathematica. Dans [Ioakimidis 93], le

logiciel Mathematica est utilisé dans la résolution des problèmes simples d‟élasticité linéaire.

Le principe de l‟approche est de conserver un paramètre sous forme symbolique dans les

formes matricielles élémentaires afin d‟effectuer l‟intégration numérique.

Les approches semi-analytiques numériques offrent un cadre favorable aux études

paramétriques utilisant des méthodes de résolution de type éléments finis. Mais au vue des

développements actuels des logiciels symboliques, ce type d‟approche ne peut s‟appliquer

qu‟à des problèmes de mécanique de taille restreinte.


27

2.3.2 Amélioration de la performance des codes éléments finis

Une autre application courante des outils de calcul algébrique consiste à effectuer des calculs

préliminaires sous formes symboliques dans le but d‟améliorer les performances numériques

dans lesquels ces résultats sont introduits. Dans [Yang 94], les expressions pour l‟élasticité en

statique sont évaluées algébriquement et le schéma d‟intégration est optimisé avant la

programmation manuelle des formes obtenues. Dans [Silvester 94], le schéma d‟intégration

numérique est également optimisé analytiquement, en utilisant le système Maple, mais du

code élément finis en FORTRAN est directement généré grâce à l‟intermédiaire d‟une

fonctionnalité de Maple. Une approche similaire est développée dans [Yagawa 90].

L‟utilisation de logiciels de calcul algébrique a montré d‟une part, qu‟il est possible d‟utiliser

la puissance et la flexibilité de tels environnements pour optimiser les expressions nécessaires

à l‟évaluation des matrices élémentaires, et d‟autre part, que le code numérique peut être

généré automatiquement à partir de l‟environnement symbolique.

2.3.3 Efficacité et flexibilité dans le développement de codes éléments finis

Une première illustration est donnée dans [Gunderson 71]. Ce document présente les

principales caractéristiques nécessaires pour développer des matrices de rigidité d‟éléments

finis à l‟aide d‟un ordinateur. Une illustration est faite sur le développement d'un élément en

plaque triangulaire troisième ordre de pliage. Cela permet de tester, à faible coût, de nouveaux

éléments pour résoudre un problème pratique donné. [Hoa 80] et [Cecchi 77] se sont basés sur

la même approche. Dans [Barbier 92], le package mathématique REDUCE est utilisé pour

produire automatiquement les matrices élémentaires de masse et de rigidité en utilisant des

polynômes d'Hermite, puis générer le code éléments finis correspondant en FORTRAN. De la

même manière dans [Korncoff 79], la production symbolique de la matrice de rigidité est

atteinte. Ici, les auteurs ont tiré profit des capacités de MACSYMA : une bibliothèque qui

donne accès à un ensemble de formes de matrices prédéfinies pour les propriétés du matériau

en élasticité linéaire. Dans [Noor 81], le potentiel d'utilisation de la manipulation symbolique

dans le développement d'éléments finis non linéaire est montré. C'est le premier travail qui a

été conçu pour traiter des problèmes non linéaires. Dans cet article, le développement

d'éléments finis non linéaires passe par trois étapes: génération des expressions algébriques

pour les coefficients de la matrice de rigidité, la génération du code correspondant en

FORTRAN pour l'évaluation numérique des coefficients et la vérification de la cohérence du

code FORTRAN généré en la comparant au tableau de coefficients donnés dans le format

MACSYMA. Dans [Cameron 97], le logiciel de calcul algébrique Maple est utilisé pour le


28

calcul des polynômes multivariés pour des modèles éléments finis. Les polynômes et leurs

dérivés sont calculés en utilisant la méthode de Hörner et des codes efficaces en C et en

FORTAN sont produits. Dans [Leff 91], un système de génération de matrice de rigidité

globale est décrit. Un fichier d'entrée pour un problème spécifique est créé pour un système

appelé SFEAS (Symbolic Finite Element Analysis System) qui génère un fichier dans le

langage symbolique mathématique REDUCE.

Toutefois, ces approches n'offrent pas forcément toute la flexibilité pour s‟attaquer aux défis

des problèmes multi-physiques et multi-échelles qui se posent dans l'industrie et la recherche.

De nouvelles approches basées sur le calcul symbolique ont été proposées afin d‟automatiser

l'élaboration de modèles éléments finis, en particulier pour la partie de discrétisation. Ces

approches consistent généralement à construire automatiquement les matrices élémentaires et

à les intégrer dans des codes plus classiques ou dans des bibliothèques d'éléments finis.

[Eyheramendy 96 et 98] propose une nouvelle manière de programmer et de développer des

modèles éléments finis en s‟appuyant sur une approche hybride symbolique/numérique pour

la résolution de problème en mécanique et sur un outil informatique de haut niveau, la

programmation orientée objet (langages Smalltalk et C++). Il présente un environnement

symbolique (FEM_Theory) pour le développement de code éléments finis qui comporte

quatre modèles (voir Figure 1). Le premier permet de représenter symboliquement les

équations aux dérivées partielles Le deuxième permet de construire et manipuler des

formulations variationnelles. Le troisième permet de construire automatiquement la forme

discrétisée du problème qui dans un quatrième modèle est intégré automatiquement dans un

code éléments finis. Cet environnement est capable de gérer tous les concepts nécessaires à la

solution de problèmes physiques : manipulations d‟équations aux dérivées partielles, formes

variationnelles, intégration par parties, formes faibles, approximations E.F... Le résultat de ces

opérations algébriques est un ensemble de données élémentaires (matrices de rigidité, de

masse, tangentes, ...) à introduire dans un code numérique classique. Le potentiel de

l‟approche est mis en évidence, d‟une part par la variété des problèmes abordés, aussi bien

dans le domaine de la mécanique linéaire (élasticité en dynamique, thermique,...), que non-

linéaire (problèmes à convection dominante, écoulement de Navier-Stokes), et d‟autre part par

le type de formulations manipulées (formulations de Galerkin, formulations de Galerkin

espace/temps continues en espace et discontinues en temps, formulations de Galerkin

stabilisées par ajout de termes de type moindres-carrés).


29

Figure 1 L‟approche FEM_Theory

[Korelc 97 et 02] a présenté une approche basée sur le logiciel Mathematica permettant de

générer automatiquement des formulations éléments finis (formes élémentaires et aspects

constitutifs). Dans un premier temps, il présente la bibliothèque SMS (Symbolic Mechanics

System) pour la dérivation automatique des formules utilisées pour l'analyse éléments finis

non-linéaire. Il présente ensuite deux bibliothèques. La première, AceGen, est un outil de

génération automatique de code. La seconde Computational Templates, combine le système

de calcul Mathematica avec une technique de Différentiation Automatique (AD). Les

caractéristiques générales d‟AceGen sont: simplification simultanée des expressions,

utilisation de techniques de différenciation automatique pour la génération automatique de

code, sélection automatique des variables intermédiaires, optimisation de procédures basées

sur l‟évaluation d‟expressions intermédiaires, lnterface d‟Utilisateur avancée, génération de

code multi-langage (FORTRAN, C/C++, Mathematica) (voir Figure 2).

Figure 2 L‟approche globale présentée par [Korelc 97, 02]

Les deux bibliothèques AceGen et Computational Templates permettent dans un contexte

multi langages de générer un code éléments finis non-linéaire à partir d‟une même description

symbolique dans des langages différents (FORTRAN, C/C++, Mathematica). L'approche

développée dans AceGen tente d‟éviter le problème habituel de la croissance des expressions

en combinant plusieurs techniques: les capacités symboliques et algébriques de Mathematica,

une technique de différentiation automatique, une génération automatique du code intégrant

Forme Forte

Lu=f

Forme

variationnelle

a(u,v)+b(v)=0

Forme

semi-discrète

Kd=f

Code

élément fini

1-SMSInitialize

2-SMSModule

3-SMSReal

4-SMSTypical

5-SMSExport

6-SMSWrite

AceGen

Code

generator

Computational

Templates

Packages

FE environnement

C Mathematica FORTRAN


30

l'optimisation simultanée d'expressions. L‟approche a été conçue pour des problèmes

complexes multi-champs. Les modules de Computational Templates sont une extension du

générateur de code AceGen. Il permet la génération de codes éléments finis multi-

environnements issue de la même description symbolique. Le temps de dérivation du code est

faible par rapport au temps de simulation typique des problèmes de taille industrielle.

[Logg et al, 07] ont développé le projet FEniCS pour l‟automatisation de la modélisation

éléments finis. FEniCS est un projet qui est constitué de plusieurs environnements (voir

Figure 3) :

- FIAT (FInite element Automatic Tabular) est un tableur d‟élément fini. Son rôle consiste

principalement à automatiser la génération des fonctions de base des éléments finis,

simplifier la spécification de nouveaux éléments et traduire les conditions sur l'espace des

fonctions en relations d‟algèbre linéaire. Il est implémenté en Python.

- FFC (FEniCS Form Compilator) est un compilateur de forme variationnelle qui automatise

l‟étape clé dans l‟implémentation de la méthode des éléments finis pour la résolution des

équations aux dérivées partielles. Il prend en entrée un ensemble d'espaces de fonctions

discrètes avec une forme multilinéaire définies sur ces espaces et produit en sortie du code

optimisé correspondant aux formes variationnelles.

Figure 3 L‟environnement FEniCS

Enfin, le code Comsol Multiphysics suit en quelque sorte une stratégie similaire, mais les

algorithmes liés à la résolution éléments finis sont figés et donc imposés pour la résolution. Il

est pourtant souvent indispensable d‟avoir accès à tous les niveaux de l‟algorithmique dans la

construction de schémas de résolution de nouveaux problèmes, notamment dans une phase de

débogage. Dans cet outil, les équations du problème peuvent être introduites soit sous forme

forte, soit sous forme faible. L‟avantage de l‟outil est qu‟il offre une infrastructure

relativement conviviale pour gérer la globalité d‟un problème depuis les équations du

problème jusqu‟à la définition de la géométrie, des actions extérieures,…

Forme

multilinéaire

Parser FFC

FIAT

Code.h


31

2.4 Une approche orientée objet intégrée symbolique/numérique

pour l’analyse éléments finis

L‟objectif de ce travail est de construire un système pour un développement rapide de

formulations éléments finis incluant le traitement de formulations mixtes complexes. Les

dérivations traditionnellement effectuées à la main sont effectuées au sein d‟un

environnement informatique de manipulations algébriques. L‟intégration entre les formes

mathématiques du problème de base (forme variationnelle) et le code calcul éléments finis est

réalisée par un environnement pour la dérivation symbolique s‟appuyant sur un concept objet

(voir Figure 4).

Figure 4 Positionnement de l‟environnement symbolique/numérique

Comme nous l‟avons montré dans le paragraphe précédent, les principales approches

développées jusqu‟à aujourd‟hui (Logg, Korelc, Eyheramendy) pour assister le

développement de codes éléments finis se basent sur des approches variationnelles pour la

description des problèmes éléments finis. Ce cadre tout à fait général est en effet aujourd‟hui

largement admis et incontournable en mécanique. Dans l‟approche proposée dans

[Zimmermann et Eyheramendy 96] et [Eyheramendy 98, 01], un environnement permettant de

gérer le problème à partir de sa forme forte sous forme symbolique est proposé.

L‟environnement permet alors de construire et manipuler la forme variationnelle pour obtenir

la forme faible. Il ne nous apparaît pas indispensable de permettre ce niveau de traitement

symbolique, à moins de proposer des outils extrêmement avancés pour couvrir le champ des

possibles des méthodes numériques développées aujourd‟hui. Cette voie reste cependant à

évaluer pour le futur. Par contre, les aspects liés à la discrétisation automatique de formes

variationnelles sont indispensables pour permettre un développement rapide de code de

simulation. Dans [Logg 07], la formulation variationnelle est introduite avant d‟être

interprétée et compilée pour générer le code correspondant. Cette description s‟appuie sur

l‟analyse tensorielle mais le formalisme ne va pas au delà de cette description, la phase de

pré-compilation de la forme variationnelle générant le code correspondant aux tenseurs et

Code

élément fini

Environnement

symbolique/numérique

Forme

variationnelle


32

opérateurs différentiels appliqués à ceux-ci. Le cadre formel défini dans ce travail ne nous

semble pourtant pas assez général. Il nous semble que pour gagner en généricité, il est

nécessaire de permettre d‟effectuer formellement les calculs des grandeurs élémentaires

quelque soit les expressions contenues dans la forme variationnelle, surtout si on souhaite

étendre l‟approche dans un contexte matériau. Enfin, il semble que dans cette approche

l‟environnement n‟est structuré que par une succession d‟outils fonctionnels et non par un

cadre représentant le modèle mathématique comme par exemple dans [Eyheramendy 96].

Dans [Korelc 97], l‟environnement développé s‟appuie sur une approche tensorielle et les

calculs formels sont réalisés dans l‟outil Mathematica. On peut noter que l‟approche s‟étend

aux formulations des modèles de matériaux. Tout comme dans [Logg 07], l‟approche est

structurée par un ensemble d‟outils fonctionnels. Ce type d‟approche ne fournit qu‟un cadre

peu structurant pour l‟approche mathématique de l‟outil développé. En effet d‟un point de vue

génie logiciel, cela constitue un handicap à long terme dans un contexte dans lequel les

évolutions en matière de modèles physiques et de méthodes numériques sont rapides. La

résolution de problématiques nouvelles basées sur des méthodes nouvelles peut entraîner une

mise en cause globale de l‟outil, et donc un redéveloppement complet ou presque.

Dans ce contexte, nous proposons une approche basée sur une description variationnelle des

problèmes de la mécanique. Le traitement des termes de la forme variationnelle doit être

suffisamment générique et donc faire appel à des techniques de calcul formel comme proposé

dans [Zimmermann & Eyheramendy 96] ou de manière plus générique dans [Korelc 02] grâce

à Mathematica. C‟est une condition qui permet d‟envisager le traitement de tout type de

modélisation physique issue de l‟industrie ou de la recherche. Enfin, contrairement aux

approches de [Logg 07] et [Korelc 97], nous allons proposer une structuration logicielle de

haut niveau d‟abstraction, appliquée au problème. Un cadre à objets doit permettre au modèle

logiciel de représenter au mieux le modèle mathématique décrivant le problème. L‟approche

est construite en deux temps. Dans un premier temps, les formes élémentaires d‟un problème

sont construites à partir de la description de l‟approche variationnelle. Dans un second temps,

ces contributions élémentaires, obtenus sous formes symboliques sont intégrées

automatiquement dans un code éléments finis.

Cette approche consiste à construire le problème discret à partir d‟un problème variationnel

quelconque en effectuant un traitement symbolique liée à la forme (voir Figure 5). La forme

variationnelle est ici discrétisée, et les matrices élémentaires sont obtenues en effectuant un


33

produit des tenseurs représentant le problème. La seconde partie de l‟approche consiste à

intégrer automatiquement les formes élémentaires dans un code de calcul comme le montre la

Figure 6. Enfin, l‟approche proposée est intégrée dans Java. Les différents temps de

l‟approche sont ainsi structurés en Java :

- Description de la formulation variationnelle.

- Discrétisation de la forme variationnelle.

- Intégration des formes élémentaires dans un environnement de simulation.

Figure 5 Une approche tensorielle pour la dérivation de modèles éléments finis

Figure 6 Intégration des formes élémentaires dans un code de calcul éléments finis

Le dernier point à souligner est le choix d‟avoir une approche logicielle complètement

homogène dans toute la phase de traitement du problème : équations du problème par

approche variationnelle, code de simulation éléments finis. En effet, nous pensons qu‟à

moyen terme, un des enjeux dans le développement des outils de calcul est l‟intégration des

aspects équationnels (descriptions de problèmes nouveaux) au sein de l‟outil de simulation et

Code EF

Algorithme de solution :

Statique

Construire le problème discret

Calcul de

Assemblage de

Calcul de

Assemblage de

Résoudre

Traitement

symbolique

Forme

variationnelle

Contributions

élémentaires

-Forme bilinéaire

-Représentation tensorielle

-Discrétisation

-Produit tensoriel

Formalisme

tensoriel

Code EF


34

de leur algorithmique. L‟outil de simulation Comsol [Littmarck 98] en est une bonne

illustration. Les outils modernes de calcul ne se contentent plus d‟intégrer des modèles

d‟équations particulières mais véritablement des méthodes de traitement de problème, tant du

point de vue équationnel que du point de vue traitement algorithmique, et ce d‟une manière

transparente pour l‟utilisateur final.

Chapitre 3 Une approche tensorielle générique pour la dérivation des modèles éléments finis

35

Chapitre 3 Une approche tensorielle

générique pour la dérivation des

modèles éléments finis

Dans ce chapitre, nous développons une approche générique dans le but de décrire la

dérivation de formulations éléments finis. Il s'agit de poser les bases d'un environnement

logiciel basé sur la dérivation de formes variationnelles ici dans un contexte éléments finis, et

permettant le développement automatisé d'outils de simulation. Nous intégrons dans

l‟approche tous les aspects nécessaires au traitement de problèmes couplés en mécanique.

Nous décrivons les différentes étapes de l‟approche. Dans le paragraphe 3.1, nous présentons

le formalisme mathématique sur lequel nous développons ce cadre formel. Dans le paragraphe

3.2, nous illustrons le formalisme mathématique sur deux exemples simples d‟équations aux

dérivées partielles : l‟élasticité linéaire et une équation non linéaire.


36

3.1 Description mathématique de l’approche de discrétisation

éléments finis

3.1.1 Formulation variationnelle

Soit

un ensemble d‟espaces des fonctions solutions et

un ensemble

d‟espaces des fonctions tests où est le nombre des espaces considérés. On considère les

formes et définies respectivement sur les produits des espaces et

et on considère le problème variationnel posé sous la forme générique

suivante:

où est le nombre de termes dans la forme variationnelle , est le ème terme de cette

forme, est l‟index de la fonction test correspondant au ème terme, est le nombre de

termes de la fonctionnelle , est le ème terme de cette fonctionnelle et est l‟index de

la fonction test dans le ème terme.

On peut toujours ramener les formes et à des formes bilinéaires et linéaire (par une

technique de linéarisation si la formulation n‟est pas linéaire) et le problème peut alors

s‟exprimer sous la forme:

où

est une forme bilinéaire,

est une forme linéaire.

Et et sont les nombres de termes respectivement des membres de droite et de gauche de

l‟équation.

3.1.2 Discrétisation éléments finis

La méthode des éléments finis appliquée à l'équation (1) consiste à remplacer par une

paire d‟espaces construits à partir de fonctions continues par morceaux.


37

On considère les approximations éléments finis classiques définies sur un élément (voir

Figure 7):

où

et

sont les fonctions d‟interpolation sur l‟élément et

et

sont les degrés de liberté correspondant respectivement à et .

Figure 7: Maillage éléments finis, élément et élément de référence

On peut donc écrire pour chacun des termes du membre de droite de la formulation

variationnelle de (1) :

Où

et

sont respectivement les vecteurs élémentaires des inconnues nodales et

des inconnues virtuelles nodales, est la matrice élémentaire correspondant à la

discrétisation du terme

, représente les interpolations des différents champs.

3.1.3 Représentation tensorielle

L‟opérateur appliqué aux fonctions d‟interpolations provient du produit de termes de la

forme variationnelle (un produit de champs auxquels sont appliqués des opérateurs

différentiels). On peut montrer que se calcule explicitement en considérant les produits

K


38

des champs (exprimés sous forme tensorielle) qui composent le terme variationnel. Il peut

alors s‟exprimer sous la forme suivante :

avec , représente le nombre de fonctions de base dans la discrétisation du

champ s, est le nombre de facteurs dans la forme , est l‟operateur différentiel

appliqué sur le champ, représente l‟ensemble des indices du terme autre qu‟un champ test

ou un champ solution, correspond aux degrés de liberté, correspond au système de

coordonnées de bases et représente l‟ensemble des indices dans la fonction discrétisée.

Cette forme est assez générale pour représenter tout type de termes de formes variationnelles

(forme multilinéaire ou forme non linéaire) dont les termes de base sont des champs

tensoriels. Cette analyse nous servira de base dans l‟implantation sous forme symbolique de la

méthode de discrétisation des termes de la forme variationnelle.

Il est important de noter que l‟approche tensorielle sur laquelle est basé le modèle est valide,

quelque soit le système de coordonnées choisi, à condition que la base de travail soit

orthonormée (voir par exemple [Garrigues 07]). Il est simplement nécessaire de donner les

expressions des opérateurs différentiels dans le système de coordonnées, une fois que celui-ci

est choisi. A ce niveau, on peut étendre l‟écriture de ces tenseurs élémentaires aux repères non

orthonormés et traiter de manière formelle les métriques (voir [Garrigues 07]).

3.2 Illustration du formalisme sur des formulations simples

3.2.1 Application à l’élasticité linéaire

On considère la forme classique d‟un problème d‟élasticité linéaire en 2D :

Avec est le vecteur déplacement, est la fonction test, est le tenseur constitutif d‟ordre 4

de la forme : , est une force et est le tenseur de

déformation défini par :

.


39

Dans ce cas trivial, les termes de la forme variationnelle sont linéaires. On pose donc

l‟inconnue et on obtient après la décomposition suivant l‟équation (1), pour et

:

On a un seul terme dans la forme, on peut donc noter que :

Donc la formulation devient:

avec :

L‟interprétation de la matrice élémentaire en appliquant l‟équation (2) devient dans ce cas :

On peut alors identifier les différents termes de la relation :

, , ,

,

,

,

;

;


40

On a donc:

où les crochets représentent les tenseurs provenant de la discrétisation de la forme

variationnelle.

En effectuant le produit doublement contracté entre les deux premiers tenseurs d‟ordre

et

suivant les indices et , on obtient le tenseur d‟ordre

qui est doublement contracté avec le tenseur d‟ordre

suivant les

indices et pour obtenir le tenseur d‟ordre 2 qu‟on le note .

La dimension de chaque tenseur est identifiée en s‟appuyant sur l‟identification des

paramètres ci-dessus. On a donc :

Pour représenter un tenseur d‟ordre supérieur à deux d‟une manière explicite, on note pour

chaque composante l‟indice correspondant. On peut donc écrire le tenseur d‟ordre 3

de dimensions sous cette forme :


41

On donne les composantes non nulles de tenseur constitutif d‟ordre de

dimensions sous cette forme :

.

Le tenseur de dimensions issu de la première contraction s‟écrit sous

cette forme :

Puis le tenseur d‟ordre 2 de dimensions s‟écrit sous cette forme :

3.2.2 Application à un problème non linéaire simple

On considère l‟équation différentielle non linéaire suivante :

Trouver fonction scalaire tel que :

La condition de Dirichlet signifie qu‟une valeur est prescrite pour l'inconnue sur un

sous-ensemble du bord du domaine et est une fonction scalaire donnée.

On pondère chacun des membres de l‟équation par la fonction test et on intègre les deux

membres sur le domaine pour obtenir le problème variationnel après intégration par parties:

Trouver fonction scalaire tel que fonction scalaire :

On défini l‟espace des fonctions test et l‟espace des fonctions solutions respectivement par :


42

et .

On peut alors écrire l‟équation (3) sous la forme suivante en utilisant le formalisme de

l‟équation (1) avec et , on obtient :

où :

Le problème est un problème non linéaire que l‟on va linéariser en utilisant une technique

classique. On écrit le développement de Taylor de au voisinage de dans la direction :

est la dérivée de évaluée en , est un produit scalaire, est l‟incrément et

est tel que

Si on néglige le reste , on définit alors la fonction comme la partie linéaire de G au

voisinage de u pour la direction :

La fonction est la linéarisation de la fonction au point et sa

dérivée directionnelle au point dans la direction .

On peut calculer la dérivée directionnelle de en dans la direction par la formule la plus

générale :

Donc finalement on obtient la linéarisation de au point :


43

On applique cela aux deux fonctionnelles :

Le calcul des dérivées directionnelles donne :

On obtient enfin les deux formes bilinéaires :

L‟équation linéarisé devient :

On pose et l‟équation s‟écrit alors :


44

Cela correspond bien à une forme bilinéaire en et pour le membre de droite de l‟équation

et une forme linéaire en pour le membre de gauche de l‟équation. On représente cette forme

bilinéaire élémentaire sous la forme :

On considère d‟abord le premier terme :

En utilisant le formalisme de (2) pour discrétiser le terme l‟équation précédente en

considérant ce premier terme :

On identifie les indices:

, , , ,

, , ,

;

;

Donc:


45

Le tenseur est un tenseur d‟ordre de dimensions

Le tenseur d‟ordre 2 de dimension

De la même manière, on considère le second terme


, , , ,

, , ,

;

;


46

Donc:

On trouve les mêmes tenseurs et

que pour

Enfin, le tenseur d‟ordre 2 de dimension :

Avec

On considère le dernier terme


, , , ,

, , ,

;

;

Donc:


47

Le tenseur de dimension est de la forme

Enfin, le tenseur d‟ordre 2 de dimension :

Finalement, on additionne les trois tenseurs pour obtenir le tenseur élémentaire

A noter que ce formalisme mathématique présenté est générique pour toute forme

variationnelle.

3.3 Analyse préliminaire sur l’implantation logicielle de la

discrétisation des formes variationnelles

L‟implantation logicielle peut être conduite selon deux voies principales. La première consiste

à construire le tenseur élémentaire de manière formelle. Cette approche a par exemple été

adoptée par [Korelc 02]. Les formes sont construites sur une approche similaire dans

Mathematica. L‟avantage de cette approche est son caractère générique. On peut


48

effectivement appliquer cela à tout type de formulation variationnelle (quelque soit le type de

termes mis en jeu). Un des inconvénients de cette approche est qu‟il est nécessaire de définir

des procédures pour obtenir toutes les formes de matrices élémentaires. De manière générale,

il faut être conscient que ce type d‟approches formelles peut conduire à l‟obtention

d‟expressions de grande taille difficiles à simplifier. La deuxième approche consiste à

élaborer des formes génériques permettant de reproduire les différents types de termes

variationnels, et de les préprogrammer. On obtient alors la forme du terme par une sorte

d‟édition de lien. Cette approche permet d‟avoir des formes optimisées des expressions pour

le calcul numérique. L‟inconvénient de cette approche est que la gestion des différents types

de termes de formes variationnelles doit être anticipéé, alors que dans la première approche

seules les expressions des opérateurs différentiels doivent être décrites et donc anticipées.

Cette approche a par exemple été adoptée par [Logg 07]. Dans notre approche, nous

choisissons de nous appuyer sur le premier type afin d‟obtenir sous forme symbolique les

formes élémentaires, et ainsi, avoir un caractère général sur cette construction. Par contre,

nous choisissons d‟implanter la démarche dans un environnement structuré afin de consolider

les approches développées.

Chapitre 4 Un environnement orienté objet pour la dérivation symbolique de formulations

éléments finis

49

Chapitre 4 Un environnement orienté

objet pour la dérivation symbolique de

formulations éléments finis

4.1 Un environnement à objets pour le calcul formel

La démarche proposée pour la dérivation de formulations éléments finis est illustrée en

Figure 8 sur l‟équation de l‟élasticité linéaire. Partant de l‟écriture d‟une formulation

variationnelle donnée, les équations de ce problème sont discrétisées. Les contributions

élémentaires du problème sont ainsi obtenues. A ce stade, le code correspondant est intégré au

sein d‟un code éléments finis classique.

L‟environnement développé intègre chacune de ces trois étapes et implémente les opérations

de discrétisation et de génération de code automatique. Nous allons donc développer un

environnement qui doit être capable de :

Représenter des formes variationnelles

Représenter des contributions élémentaires et de les gérer dans le contexte de la

formulation

Réaliser des calculs formels pour obtenir ces contributions élémentaires


éléments finis

50

Générer du code de calcul dans un langage donné, pour une structure de codes

éléments finis donné, à partir de ces contributions élémentaires.

Figure 8 Approche variationnelle pour modèles éléments finis

Cette approche sous-entend le développement d‟un moteur de calcul formel, y compris

matriciel, afin d‟obtenir sous forme symbolique toutes les contributions élémentaires. Un des

points clés de l‟approche que nous proposons est de n‟utiliser qu‟un seul et même paradigme

de programmation, orienté objets en Java, de façon à tendre vers une intégration complète et

naturelle des outils de simulation depuis les équations du problème jusqu‟à la simulation d‟un

problème particulier. L‟avantage d‟une telle démarche est qu‟elle permet d‟envisager à terme

d‟intégrer le développement d‟algorithme de résolution numérique de l‟outil de simulation

(code éléments finis) sous forme symbolique, et donc d‟ouvrir un espace de liberté dans

l‟élaboration du code de simulation final.

Formulation variationnelle

,

+ Conditions limites + Conditions aux bords

Contributions élémentaires

Code élément fini

Discrétisation automatique :

- En temps et en espace

- Analyse tensorielle

- …

Génération automatique :

Génération automatique d‟un

élément dans un code existant


éléments finis

51

La structure de l‟environnement est présentée dans ce chapitre. Dans le paragraphe 4.2, les

objets fondamentaux nécessaires pour la représentation symbolique d‟une forme

variationnelle sont présentés. La hiérarchie des classes issue de cette étude y est présentée.

Les objets nécessaires pour l‟automatisation de la discrétisation sont présentés dans le

paragraphe 0. L‟implantation est réalisée Java.

4.2 Analyse à objet pour la représentation symbolique d’une

forme variationnelle

L‟objectif premier dans l‟environnement symbolique est de construire un environnement

capable de représenter les équations d‟une formulation variationnelle. Cela fait intervenir des

champs scalaires, vectoriels et tensoriels comme l‟illustre la Figure 9 sur le cas de l‟élasticité

linéaire en petites déformations/petits déplacements : u et v sont des champs vectoriels de

déplacement respectivement solution et test, est un tenseur d‟ordre 2 (opérateur

différentiel appliqué à u), est un tenseur d‟ordre 2 (opérateur différentiel

appliqué à v), C est le tenseur constitutif d‟ordre 4, f est un vecteur de même dimension que u

et v. Les produits et sommes de termes sont intégrés sur un domaine, les intégrales sont

sommées et forment les membres de droite et de gauche des équations.

+ Conditions limites + Conditions aux bords

Figure 9 Exemple de forme variationnelle sur l‟élasticité linéaire

L‟objet le plus naturel à définir est le champ. On le définit par son nom. C‟est une grandeur

tensorielle qui peut être défini par sa dimension, son ordre, et qui peut connaître des données

supplémentaires telles que les fonctions d‟interpolations,… Deux exemples sont donnés dans

la Figure 10. Le champ est un tenseur d‟ordre 1 et dimension 2 (vecteur de dimension 2)

qui sera interpolé par les fonctions sur un élément de référence, comme classique en

élément finis (voir [Hughes 87] par exemple).

Le tenseur constitutif , où et sont les coefficients du

Lamé, est un tenseur d‟ordre 4. Il n‟est pas un champ solution ou test, et il n‟est pas interpolé.


éléments finis

52

Champ :

nom: u

dimension : [2] // Tenseur d‟ordre 1,

dimension de l‟espace: 2

type de champ: solution

fonctions d‟interpolation : [N1, N2 N3, N4]

Champ:

nom: C

dimension : [2,2,2,2] // Tenseur d‟ordre 4

dimension de l‟espace : 2

type de champ: -

fonctions d‟interpolation: -

Figure 10 Exemples d‟objets champ

On peut maintenant définir le tenseur des déformations simplement en appliquant au

champ l‟opérateur différentiel (.) comme le montre la Figure 11. L‟objet opérateur

différentiel est ici appliqué à un champ. La principale opération qu‟il réalise est le calcul de

l‟expression correspondante à l‟opérateur appliqué au champ.

Opérateur différentiel :

nom :

champ :

Figure 11 L‟opérateur différentiel

Les composantes de l‟opérateur différentiel sont des dérivées partielles. Ces dérivées

partielles dépendent du champ et de l‟indice de la variable de dérivation, comme on le montre

dans l‟objet dérivée partielle dans la Figure 12.

Dérivée partielle :

champ :

indice de dérivation : 1

Figure 12 L‟objet dérivée partielle

Pour exprimer les composantes de cet objet, nous avons besoin de l‟objet tenseur. Cet objet

donne la représentation tensorielle d‟un champ ou d‟un opérateur différentiel appliqué à un

champ.

Dans la Figure 13, nous donnons la représentation tensorielle de l‟opérateur différentiel

en le définissant comme un tenseur d‟ordre 2, de dimension 2 et ses composantes sont

calculées dans l‟objet dérivée partielle.


éléments finis

53

Tenseur :

dimension : 2

ordre : 2

composantes :

Figure 13 L‟objet tenseur

On introduit un nouveau type d‟opérateurs permettant de représenter les opérations entre

champs telles que les opérateurs binaires : la somme, le produit scalaire, le produit tensoriel,

le produit contracté,…. et les opérateurs unaires : intégrale,…Ces opérateurs permettent

d‟effectuer les différentes opérations algébriques entre les champs pour représenter la forme

variationnelle.

L‟objet produit doublement contracté effectue une double contraction entre deux tenseurs

d‟ordre supérieur ou égal à deux. Dans la Figure 14, on effectue une double contraction entre

le tenseur d‟ordre 2 et le tenseur d‟ordre 4 .

Contraction double :

membre de gauche :

membre de droite :

Figure 14 L‟objet produit doublement contracté

L‟objet produit scalaire effectue un produit scalaire entre deux vecteurs ou entre deux

scalaires ou entre un vecteur et un scalaire. Dans la Figure 15, on effectue un produit scalaire

entre le vecteur et le vecteur .

Produit scalaire :

membre de gauche :

membre de droite :

Figure 15 L‟objet produit scalaire

Une ou plusieurs opérations effectuées génèrent une expression. L‟objet expression représente

un terme composé d‟opérateurs et de champs. Dans la Figure 16, on effectue deux doubles

contractions pour obtenir une expression.


éléments finis

54

Expression :

opérateurs : et :

champs : et

Figure 16 L‟objet expression

L‟objet intégral effectue l‟intégrale d‟une expression sur un domaine défini. Dans la Figure

17, on représente l‟intégrale de l‟expression sur le domaine . Cette intégrale

est aussi une expression.

Intégrale :

domaine :

expression :

Figure 17 L‟objet intégrale

L‟objet somme effectue une somme entre deux expression. Dans la Figure 18, on effectue la

somme entre deux intégrales.

Somme :

gauche de l‟opérateur :

droite de l‟opérateur :

Figure 18 L‟objet somme

L‟objet équation est crée d‟une manière évidente. Cet objet contient une expression dans

chaque membre (droite et gauche). Il est illustré dans la Figure 19.

Equation :

membre de gauche :

membre de droite :

Figure 19 L‟objet équation

La Figure 20 met en évidence l‟implantation de la structure d‟objet pour représenter

l‟équation de l‟élasticité linéaire. Les objets opérateurs sont représentés par des rectangles, les

champs par des ellipses. C‟est une structure arborescente qui permet de représenter tout type

d‟équations.


éléments finis

55

Figure 20 Hiérarchie d‟opérateur.

La quasi-totalité des informations nécessaires au développement du code correspondant à la

formulation variationnelle sont contenues dans l‟objet formulation. Cet objet regroupe tous les

éléments d‟une formulation variationnelle. L‟objet formulation contient : les équations du

problème, les contributions élémentaires (issues de la discrétisation) et les caractéristiques du

matériau utilisé (voir Figure 21) et bien sûr connaît les champs inconnus et champs tests de la

formulation.

Formulation : Elasticité linéaire

-équations :

-contributions élémentaires :

-caractéristiques du matériau : ,

Figure 21 L‟objet formulation

Le rôle principal de l‟objet formulation est d‟intégrer ses contributions (contributions

élémentaires du problème éléments finis) dans le code de simulation numérique.

L‟analyse effectuée dans le paragraphe précédent a conduit à définir un certain nombre

d‟objets. Une analyse supplémentaire a conduit à définir les classes correspondantes sont

hiérarchisées comme le montre la Figure 22. L‟organisation hiérarchique est conçue de façon

à permettre la construction des structures arborescentes pour représenter les formes

=

+

ʃ(.)dx

:

:

ε(.)

u

C

ε(.)

v

ʃ(.)dx

.

f v

0

: Opérateur

: Champ


éléments finis

56

variationnelles et autres expressions nécessaires. Nous avons adopté un formalisme UML

(Unified Modeling Language) pour cette figure.

Figure 22 Hiérarchie partielle des classes

La première classe, base de toutes les autres classes dans la Figure 22, est la classe Operator.

Cette classe est à la base de la structuration des expressions et gère les différentes opérations.

Son existence permet de créer des structures d‟expressions arborescentes.

Operator

-name : String

-leftConnected :Operator

-rightConntected :Operator

+sum()

+minus()

+product()

+substitute()

+…

Figure 23 Classe Operator

Cette classe étant la super-classe des autres classes. Cela permet de typer chaque nœud de

l‟arborescence. Les attributs de la classe sont : le nom de l‟opérateur qui est une chaîne de


éléments finis

57

caractères et deux attributs leftConnected et rightConnected qui permettent de connaître les

objets auxquels l‟opérateur est connecté et aussi de savoir s‟il est connecté sur un lien gauche

ou droit (voir Figure 23).

La classe Operator est une classe abstraite. Les méthodes sont implémentées dans les classes

dérivées. La classe UnaryOperator (voir Figure 24) sous-classe de la classe Operator

permet de représenter les opérateurs unaires. Cette classe est aussi une classe abstraite.

UnaryOperator

-operator : Operator

+putConnectedOperator()

+…

Figure 24 Classe UnaryOperator

La classe UnaryOperator permet de représenter les objets tels que les opérateurs

différentiels, les intégrales, les dérivées partielles… Afin de pouvoir intégrer les champs dans

les expressions, nous avons fait le choix de faire hériter la classe Field de la classe

UnaryOperator. La classe Field est donnée dans la Figure 25. Les attributs de la classe sont

spécifiés : le type du champ (scalaire, vecteur ou tenseur) qui est un tableau d‟entier décrivant

l‟ordre et les dimensions du tenseur, la dimension de l‟espace vectoriel dans lequel on

effectue les opérations, la nature du champ (solution ou non) et les fonctions d‟interpolations

qui forment un tableau de type Operator.

Field

-type: int []

-spaceDimension: int

-solutionOrTest: boolean

-interpolationFunction: Field[]

+substitute()

+giveString()

+…

Figure 25 Classe Field

Une instance de la classe Field est capable de s‟analyser, c‟est-à-dire de définir le type de son

champ (scalaire, vecteur ou tenseur) en fonction de son ordre et de ses dimensions définies

dans l‟attribut type. Des tâches de manipulations algébriques lui sont assignées : substitution,

addition,…


éléments finis

58

Les deux classes PartialDerivative et DifferentialOperator permettent de représenter

respectivement les dérivées partielles appliquées à une fonction et les opérateurs différentiels

appliqués à un champ. La structure et le comportement de ces classes sont les mêmes que

ceux de la classe UnaryOperator. Pour la première, nous avons deux attributs

supplémentaires qui correspondent : l‟un à l‟indice de dérivation et l‟autre à l‟ordre de

dérivabilité en temps. On récupère la chaine de caractère qui correspond au terme (méthode

giveString()) (voir Figure 26).

PartialDerivative

-indice: int

-timeOrderDerivative:int

+giveString()

+substitute()

+expandAll()

+…

Figure 26 Classe PartialDerivative

Pour la classe DifferentialOperator, un seul attribut supplémentaire est nécessaire, son nom.

Une instance de la classe DifferentialOperator peut déterminer l‟ordre et la dimension

correspondant au champ auquel il est appliqué, à partir des caractéristiques du champ (ordre et

dimension du champ, voir Figure 27).

DifferentialOperator

-name: String

+sum()

+expandAll()

+giveString()

+…

Figure 27 Classe DifferentialOperator

D‟un point de vue manipulations algébriques, pour un opérateur différentiel appliqué à une

expression, si celui-ci est linéaire, la propriété de linéarité est implantée dans la méthode

expandAll().

Pour définir une expression complexe, nous définissons la classe Expression. Il s‟agit de

décrire et représenter tout type d‟expression dans la forme. Cette classe dérive de la classe

UnaryOperator. Elle hérite donc de l‟attribut unique de la classe UnaryOperator qui est de

type Operator. Cet attribut constitue donc la racine de l‟arborescence représentant

l‟expression.


éléments finis

59

Figure 28 Classe Expression

La classe BinaryOperator (voir Figure 29) hérite de la classe Operator. Elle permet de

représenter les opérateurs binaires.

BnaryOperator

-left : Operator

-right : Operator


+…

Figure 29 Classe BinaryOperator

Un opérateur binaire possède deux opérandes a priori non commutatifs ; un à droite et un à

gauche. Cette classe est aussi une classe abstraite. Cette classe permet d‟effectuer les

opérations telles que : la somme, le produit scalaire,… Elle possède deux attributs de types

Operateur: left et right qui correspondent aux deux membres de l‟opérateur. On peut note

que les opérateurs binaires ne sont pas forcément commutatifs.

La classe Sum dérive de la classe BinaryOperator (voir Figure 30) représente une somme,

quelque soit le type de termes sommés (scalaire, vecteur ou tenseur). Cette classe a le même

comportement de la superclasse. Le seul attribut supplémentaire nécessaire est un attribut

binaire qui donne le signe liant les termes de gauche et de droite de la somme. La

distributivité et l‟associativité sont implantées dans la méthode expandAll() pour effectuer

l‟opération d‟expansion de l‟expression.

Sum

-sign: boolean

+distribute()

+expandAll()

+giveString()

+…

Figure 30 Classe Sum

Expression

++distribute()

+expandAll()

+giveString()

+…


éléments finis

60

La classe Product est la classe générique qui va permettre de définir tous les produits. A

l‟opposé des sommes, les produits sont différenciés suivant leur type (scalaire, produit

contractés,…)

Product

-left : Operator

-right : Operator


+…

Figure 31 Classe Product

La classe Product est donc déclarée comme classe abstraite (voir Figure 31), et les méthodes

sont implémentées dans les trois sous-classes ScalarProduct, SimpleContraction et

DoubleContraction (voir Figure 32). La classe ScalarProduct gère le produit scalaire entre

les différents types de champs tels que les scalaires, les vecteurs ou les tenseurs. Ce produit

est représenté par un espace « ». Une instance de la classe SimpleContraction gère le

produit simplement contracté deux vecteurs, deux tenseurs ou un vecteur et un tenseur. Ce

produit est représenté par un point « . ».

Figure 32 Les classes ScalarProduct et SimpleContraction

La classe DoubleContraction gère le produit doublement contracté entre deux tenseurs. Ce

produit est représenté par deux point « : » (voir Figure 33)

DoubleContraction

+distribute()

+expandAll()

+giveString()

+…

Figure 33 Classe DoubleContraction

ScalarProduct

+distribute()

+expandAll()

+giveString()

+…

SimpleContraction

+distribute()

+expandAll()

+giveString()

+…


éléments finis

61

La classe Fraction dérive de la classe BinaryOperator. Elle représente une fraction

rationnelle, dont les instances peuvent intervenir dans des expressions. Ses méthodes

permettent de gérer les manipulations algébriques au sein d‟expression complexes.

Fraction

+addition()

+division()

+multiplication()

+giveString()

+…

Figure 34 Classe Fraction

La classe Integral donnée dans la Figure 35 dérive de la classe UnaryOperator. Elle hérite

l‟attribut Operator qui est a priori d‟une expression (la racine d‟un objet Expression est

forcément un objet Operator), et possède un deuxième attribut qui décrit le domaine (de type

chaine de caractère). L‟intégrale implémente la propriété de linéarité de l‟intégrale pour

effectuer une opération d‟expansion sur une expression intégrale (méthode expandAll()).

Integral

-domain: String

+expandAll()

+substitution()

+giveString()

+…

Figure 35 Classe Integral

Une équation intégrale est une instance de la classe Equation. Cette classe dérive de la classe

BinaryOperator. Elle exprime l‟égalité entre les intégrales du membre de droite et de gauche

comme le montre la Figure 36.

Equation

+expandAll()

+substitution()

+giveString()

+…

Figure 36 Classe Equation


éléments finis

62

La classe Formulation représente la formulation variationnelle du problème. Ses attributs

sont : le nom de la formulation, les équations du problème, les contributions élémentaires du

membre de droite et de gauche et les caractéristiques du matériau (voir Figure 37).

Formulation

-name : String

-equations : Equation[]

-leftSymbolicTensor :Tensor[]

-rightSymbolicTensor :Tensor[]

-material : femjava.material

+giveElementalBlockMatrixString()

+giveFullMatrixString()

+…

Figure 37 Classe Formulation

Toutes les classes fondamentales nécessaires à la génération de formes variationnelles ont été

définies dans cette partie.

Considérons l‟exemple d‟une forme faible d‟un problème d‟élasticité donné dans la Figure 9

pour représenter une forme variationnelle. La structure arborescente est donnée dans la

Figure 38. On considère le membre de gauche. C‟est une expression qui est une intégrale

contenant une expression qui est un produit scalaire de deux champs. Le membre de droite est

une expression qui est une somme de deux intégrales. La première contient une expression qui

est un produit. Ce produit lui-même contient une expression un opérateur différentiel etc.…

Le problème est qu‟une même information peut être demandée à différents niveaux de

l‟arborescence pour effectuer différentes tâches. Il est nécessaire premièrement de ne stocker

l‟information qu‟une seule fois dans la structure, afin d‟assurer facilement la cohérence des

informations, et deuxièmement, d‟être sûr d‟y accéder lorsque nécessaire.

Cette structure permet de gérer le schéma de circulation naturelle des informations vers le bas

de l‟arborescence (lien de l‟objet avec ses attributs). Ainsi, un lien vers le supérieur

hiérarchique dans l‟arborescence rend possible la remontée des messages. Cela est implanté

dans la hiérarchie des opérateurs (voir Figure 22) en s‟appuyant sur les attributs

leftConnected et rightConnected de la classe abstraite Operator Ainsi un opérateur connaît

les opérateurs auxquels il est rattaché, qu‟il se situe dans le lien droit ou gauche de

l‟opérateur.


éléments finis

63

Figure 38 Structure arborescente pour une instance d‟Equation

Une classe supplémentaire est nécessaire pour représenter les tenseurs, matrices élémentaires

et effectuer les opérations décrites au Chapitre 3 pour construire les matrices élémentaires.

La classe Tensor permet de représenter un tenseur. Le premier attribut de cette classe est la

dimension du tenseur qui est un tableau d‟entier et qui sous-entend l‟ordre du tenseur. Le

deuxième représente les composantes. Ces composantes sont des expressions mathématiques ;

le type de l‟attribut correspondant est donc Operator. Dans cette classe, les opérations

tensorielles classiques sont implémentées : la somme (méthode sum()), le produit tensoriel

(méthode tensorProduct()), produit simplement contracté, produit doublement contracté,

transposée (méthode transpose()), inverse -si inversible- (méthode inverted()), la trace

(méthode trace()).... Le diagramme UML succinct de la classe est donné Figure 39.

Tensor

-dimension :int[]

-components :Operator[]

+tensorProduct()

+sum()

+transpose()

+inverted()

+trace()

+…

Figure 39 Classe Tensor


éléments finis

64

On peut noter qu‟une part importante du comportement de la classe consiste à gérer toutes les

opérations entre tenseur, tels que les produits entre les tenseurs. A titre d‟exemple, on

considère les produits contractés (simple, double,…). L‟ensemble des produits contractés

s‟effectue grâce à une seule méthode, la méthode tensorContract(Tensor, int[], int[]). Cette

méthode effectue le produit contracté entre le tenseur receveur et celui passé en argument, et

prend en compte les indices passés en argument pour la contraction. Les deux tableaux

d‟entiers désignent le numéro des indices à contracter respectivement pour le receveur et le

tenseur passé en argument. L‟ordre et les dimensions associées du tenseur issu d‟une

contraction est calculé à l‟aide de la méthode giveTensorContractDimension(Tensor, int[],

int[]) qui renvoie un tableau d‟entiers dont la taille est l‟ordre du tenseur, et dont les éléments

sont les dimensions repectives attachées à chaque dimension. L‟ordre de ce tenseur doit être

égal à la somme de l‟ordre de deux tenseurs contractés moins le nombre des indices contractés

(voir Figure 41).

public int[] giveTensorContractDimension(Tensor tensor,

int[] tab1, int[] tab2) {

int [] tab = null;

tab = new int[this.getOrder() + tensor.getOrder() -

tab1.length - tab2.length];

return tab;

}

Figure 40 La méthode giveTensorContractDimension(Tensor, int[], int[])

Afin d‟illustrer l‟appel de la méthode tensorContract(Tensor, int[], int[]), on considère et

deux tenseurs d‟ordre 2. Le produit simplement contracté entre et est réalisé par la

commande A.tensorContract(B, int[]{2}, int[]{1}) correspondant à l‟exemple Figure 41.

Figure 41 Produit contracté simple

Il est important de rappeler qu‟une des opérations importantes réalisée par le champ, instance

de Field, (opérateurs différentiels appliqués aux champs) est de construire sa représentation

tensorielle. Les composantes du terme considéré sont données à l‟aide de la méthode

buildTensor(). L‟ordre du tenseur obtenu est lié à la nature de l‟opérateur appliqué au champ.


éléments finis

65

Afin d‟illustrer l‟opération, on considère l‟opérateur tenseur de déformation appliqué au

champ , tenseur d‟ordre 1 et de dimension 2. Dans la Figure 42, on donne les composantes

du tenseur de dimension construit à l‟aide de la méthode buildTensor().

1,1=[du1_dx1];

1,2=[0.5*(du1_dx2 + du2_dx1)];

2,1=[0.5*(du2_dx1 + du1_dx2)];

2,2=[du2_dx2] ;

Figure 42 les composante de

4.3 Discrétisation de la forme variationnelle et formes

élémentaires

Dans ce paragraphe, nous décrivons la modélisation à objets de l‟approche mathématique de

discrétisation d‟une forme variationnelle présentée dans le Chapitre 3. Cela revient à décrire

la génération automatique de tenseur élémentaire . Nous présentons les objets nécessaires

à la construction de cette structure et également ceux nécessaires pour effectuer toutes les

opérations liées à la discrétisation.

4.3.1 Discrétisation symbolique

La première étape de la discrétisation automatique consiste à discrétiser chacun des champs

test et solution (ou incrément de solution pour les problèmes non linéaires) des termes

intégraux générés à partir des champs et opérateurs différentiels appliqués. La discrétisation

se base sur l‟approche tensorielle présentée dans le paragraphe 3.1.3 du Chapitre 3.

L‟opération de discrétisation d‟un champ auquel est appliqué un opérateur différentiel

consiste à crééer le tenseur qui représente sa discrétisation symbolique. Afin d‟illustrer

l‟approche, nous nous appuyons sur l‟exemple d‟un vecteur de dimension 2 que l‟on

interpole en 4 nœuds : fonctions d‟interpolations . On peut alors écrire

l‟interpolation du champ u sur l‟élément e de la manière suivante :

, où

les (pour et ) sont les degrés de liberté correspondant au champ u sur

un élément e.


éléments finis

66

On pose :

, vecteur des inconnues nodales sur un élément de dimension 8. On peut

alors écrire l‟interpolation du vecteur u sur l‟élément e :

On considère maintenant le tenseur des déformations . Son expression est un tenseur de

d‟ordre 2 et dimension

. On peut alors écrire la

forme discrétisée correspondant à où est un tenseur d‟ordre 3 et de

dimension . Sa notation indicielle est où est l‟indice lié à la discrétisation.

Les composantes du tenseur sont alors :

Les opérations qui correspondent à la construction de la matrice sont prise en charge par la

classe FEDiscretization qui possède un seul attribut : le champ à traiter (Figure 43).

FEDiscretization

-field

+substitute()

+giveDiscretizationTensor()

+…

Figure 43 Classe FEDiscretization

D‟un point de vue pratique, le principe de construction de la matrice est simple. En

reprenant l‟exemple du tenseur de déformations en 2D, on construit d‟abord un tenseur

d‟ordre 2+1. Pour tous les p, on affecte les coefficients . On effectue enfin une

substitution des composantes et

par les fonctions d‟interpolations correspondant

aux degrés de liberté (méthode substitute()). On obtient alors le tenseur d‟ordre par la

méthode giveDiscretizationTensor() présentée Figure 44,. Dans cette méthode (voir Figure

44), on reprend pas à pas le principe d‟élaboration de décrit ci-dessus :


éléments finis

67

- Partie A : construction de et affectation des

- Partie B : substitution des et

.

On peut ainsi traiter tous les champs tests et solutions (ou incréments de solution pour les

problèmes non linéaires) de chacun des termes intégraux d‟une formulation variationnelle.

public Tensor giveDiscretizationTensor (Field F) {

Operator[] interpFunctions=F.giveInterpolationFunctions();

//----------------------------A-----------------------------

Tensor t = new Tensor

(this.giveFETensorDicretizationDimension());

Field f=(Field) differentialOperator.giveField();

Tensor tensor=differentialOperator.getTensor();

Operator op = null;

Operator[] dofsX = ((Field) differentialOperator.

giveField()).buildBaseTensor().getComponents();

Field[] dofs = new Field[dofsX.length];

for (int ll = 0; ll < dofs.length; ll++){

dofs[ll] = (Field) dofsX[ll];

int[][] positions = tensor.givePosition();

int[] tab = null;

for (int i=0;i<tensor.giveNumberofComponents();i++) {

op = tensor.getComponent(positions[i]);

tab = new int[t.getDimension().length];

for (int k = 0; k < positions[i].length; k++)

tab[k] = positions[i][k];

}

//----------------------------B-----------------------------

int number = 1;

for (int l = 0; l < numberOfNodes; l++)

for (int j = 0; j < dofs.length; j++) {

tab[tab.length - 1] = number;

Operator copy = op.deepCopy();

Operator[] newField=new Field[dofs.length];

for (int k = 0; k < dofs.length; k++){

newField[k] = new Field("0", null);

newField[j] = interpFunctions[l];

copy=copy.substitute(dofs, newField);

t.setComponent(tab, copy);

number++;

}

}

}

return t;

}

Figure 44 La méthode giveDiscretizationTensor(Field)


éléments finis

68

L‟ordre et la dimension du tenseur sont calculés à l‟aide de la méthode

giveFEDiscretizationTensorDimension(). Cette méthode renvoie un tableau d‟entiers (voir

Figure 45).

public int[] giveFETensorDicretizationDimension() {

Tensor tensor = this.getTensor();

int[] dim = null;

dim = new int[tensor.getDimension().length + 1];

for(int i = 0 ; i < tensor.getDimension().length; i++)

dim[i] = tensor.getDimension()[i];

dim[tensor.getDimension().length] = this

.giveDiscretizationIndice();

return dim;

}

Figure 45 La méthode giveFEDiscretizationTensorDimension()

4.3.2 Génération des contributions élémentaires

La deuxième étape de la démarche du calcul de la matrice élémentaire consiste à effectuer les

produits entre tenseurs issus de la discrétisation de chacun des termes du terme intégral.

Pour effectuer ces opérations, on s‟appuie sur les sous classes de la classe abstraite Product

qui gèrent tous les types de produits entre tenseurs tels que : SimpleContraction

DoubleContraction, ScalarProduct, TensorProduct,… Dans la Figure 46, nous résumons

l‟approche globale pour obtenir la matrice élémentaire complète sur l‟exemple de l‟élasticité

linéaire. On peut donc, pour chaque terme intégral de la formulation, obtenir automatiquement

ses contributions élémentaires en effectuant les produits entre tenseurs (dans cet exemple

produits doublement contractés).


éléments finis

69

Figure 46 Obtention des contributions élémentaires

Cela est illustré sur l‟opération . Le tenseur issu de la contraction double entre

le tenseur d‟ordre 3 et le tenseur d‟ordre 4 est un tenseur d‟ordre

3 ( ) et de dimension . Cette opération est illustrée sur la Figure

47.

Figure 47 Premier produit contracté double

Le tenseur issu de la contraction double entre le tenseur d‟ordre 3 et le tenseur d‟ordre 3

est un tenseur d‟ordre 2 ( ) et de dimension .

Cette opération est illustrée sur la Figure 48.

Integral

(forme bilinéaire)

Tensor

Première étape :

représentation tensorielle

FEDiscretization

Deuxième étape :

Discrétisation

(contraction double)

Troisième étape :

Contraction

(tenseur élémentaire)

ordre : 2 ; dimension :

i, j, k, l=2 et p=8


éléments finis

70

Figure 48 Second produit contracté double

La contraction entre deux tenseurs s‟effectue à l‟aide de la méthode générique

tensorContract(Tensor, int[], int[]) décrite dans le paragraphe précédent. Les méthodes

correspondant aux produits contractés doubleContraction(Tensor),

simpleContraction(Tensor)),ect… s‟appuient sur cette méthode générale pour lequel on

définit les indices concernés par l‟opération (voir exemple de doubleContraction(Tensor T)

Figure 49).

public Tensor doubleContraction(Tensor T) {

int[] dim1 = this.giveNotLeftDiscretizationIndTab(2);

int[] dim2 = T.giveNotRightDiscretizationIndTab(2);

Tensor answer = new Tensor

(this.giveTensorContractDimension(T, dim1, dim2));

int[] tab = this.giveDiscretizationIndTab(T, this

.giveTensorContractDimension(T, dim1, dim2));

answer = this.tensorContract(T, dim1, dim2);

return answer;

}

Figure 49 La méthode doubleContraction(Tensor)

Dans la Figure 50, on donne une partie des composantes du tenseur . Ces composantes

sont données à partir de la méthode giveDiscretizationTensor().

1,1,1=[dN1_dx1*(lambda + 2.0*mu)];

1,1,2=[0.5*(dN1_dx2)*mu + 0.5*(dN1_dx2)*mu];

1,2,1=[0.5*(dN1_dx2)*mu + 0.5*(dN1_dx2)*mu];

1,2,2=[dN1_dx1* lambda];

2,1,1=[dN1_dx2* lambda.0];

2,1,2=[0.5*(dN1_dx1)*mu + 0.5*(dN1_dx1)*mu];

7,2,1=[0.5*(dN4_dx2)* mu + 0.5*(dN4_dx2)* mu];



8,1,2=[0.5*(dN4_dx1)* mu + 0.5*(dN4_dx1)* mu];

8,2,1=[0.5*(dN4_dx1)* mu + 0.5*(dN4_dx1)* mu];

8,2,2=[dN4_dx2*( lambda + 2.0* mu)]]

Figure 50 Extrait des composantes de la matrice

k, l=2 et p, q=8


éléments finis

71

Dans la Figure 51, on donne quelques composantes de la matrice . La matrice

est la contribution élémentaire du terme

.

[1,1=[dN1_dx1*( lambda + 2.0* mu)*dN1_dx1 + (0.5*(dN1_dx2)* mu

+ 0.5*(dN1_dx2)* mu)*0.5*(dN1_dx2) + (0.5*(dN1_dx2)* mu +

0.5*(dN1_dx2)* mu)*0.5*(dN1_dx2)];

1,2=[(0.5*(dN1_dx2)* mu + 0.5*(dN1_dx2)* mu)*0.5*(dN1_dx1) +

(0.5*(dN1_dx2)* mu + 0.5*(dN1_dx2)* mu)*0.5*(dN1_dx1) +

dN1_dx1* lambda *dN1_dx2];

8,7=[dN4_dx2* lambda *dN4_dx1 + (0.5*(dN4_dx1)* mu +

0.5*(dN4_dx1)* mu)*0.5*(dN4_dx2) + (0.5*(dN4_dx1)* mu +

0.5*(dN4_dx1)* mu)*0.5*(dN4_dx2)];

8,8=[(0.5*(dN4_dx1)* mu + 0.5*(dN4_dx1)* mu)*0.5*(dN4_dx1) +


dN4_dx2*( lambda + 2.0* mu)*dN4_dx2]]

Figure 51 une partie des composantes du tenseur

La quasi-totalité des informations nécessaires au développement du code correspondant à la

formulation variationnelle sont contenues dans classe Formulation. Cette contient le

problème continue (un tableau d‟équations) et le problème discret (un tableau de tenseurs) de

la formulation variationnelle. Elle contient également les caractéristiques du matériau utilisé

(voir Figure 52).

Formulation

-name : String ;

-equations : Equation[]

-leftSymbolicTensor : Tensor[] ;

-rightSymbolicTensor : Tensor[] ;

-material : Material.femjava ;

+giveInitializeMethod()

+giveFormulationName()

+giveElementalMatrices()

+…

Figure 52 La classe Formulation

Les principales méthodes de cette classe ont pour rôle de générer les entêtes des méthodes

indispensables pour la génération de code.

L‟environnement proposé permet de représenter les formes continues et les contributions

élémentaires au problème éléments finis, le tout dans un cadre formel (possibilité d‟effectuer


éléments finis

72

des traitements mathématiques). La représentation de la forme continue du problème est basée

sur une approche par champs (scalaire, vecteur, tenseur) auxquels sont appliqués des

opérateurs différentiels, les différents termes d‟une formulation variationnelle étant

représentés : sommes de produits de termes, termes intégraux,… De manière similaire, les

contributions élémentaires (matrices) sont représentées dans l‟environnement informatique.

L‟outil de discrétisation automatique qui permet de construire les contributions élémentaires

de manière formelle s‟appuie sur la description tensorielle de la forme variationnelle.

L‟environnement logiciel est structuré sur la base de cette description mathématique. Des

outils permettant d‟obtenir automatiquement les contributions élémentaires d‟un problème

donné et de les intégrer dans un code de calcul existant complètent l‟approche. Dans la suite,

on complète la procédure en intégrant ces contributions élémentaires dans un code éléments

finis.

Chapitre 5 Génération automatique de code éléments finis

73

Chapitre 5 Génération automatique de

code éléments finis

5.1 Génération automatique de code et code éléments finis

Dans les chapitres précédents nous avons développé une méthodologie et des outils pour

automatiquement générer sous forme symbolique les contributions élémentaires

correspondant à un problème variationnel donné. Sur le principe d‟un outil de pré-

compilation, il est alors très naturel de générer un code source dans un langage informatique

donné, permettant de calculer numériquement chacune de ces contributions élémentaires. Sur

le principe, cela revient à générer dans un premier temps les formules correspondant aux

composantes des contributions élémentaires, et d‟autre part, à générer les en-têtes et autres

éléments permettant l‟initialisation des différentes variables. Mais il manque cependant un

élément essentiel dans l‟approche. En effet, considérer le calcul des contributions

élémentaires ne suffit pas à gérer un calcul éléments finis global. Cela ne comprend pas par

exemple : l‟assemblage des contributions élémentaires calculées dans un système linéaire,

l‟intégration de ces contributions dans des algorithmes de résolution plus globaux (schémas

d‟intégration en temps, algorithme de résolution non linéaire,…),… Deux pistes peuvent alors

être suivies. La première, la plus simple et relativement naturelle dans notre contexte, consiste

à s‟appuyer sur un code éléments finis existant et y intégrer le code source correspondant au


74

calcul des contributions élémentaires. Dans le contexte de ce travail, il faut alors adjoindre au

code généré les éléments permettant de les intégrer dans la structure du code. La seconde piste

consiste à intégrer dans le cadre de réflexion général toute l‟approche numérique d‟un

problème variationnel, au delà de la discrétisation en espace (schéma d‟intégration en temps,

algorithme de traitement de non linéarités,…) et toute l‟algorithmique de base permettant de

gérer un calcul donné (assemblage, mise à jour des données, résolution de systèmes

linéaire,…). Cette seconde alternative dépasse le cadre de ce travail de thèse mais représente

une perspective incontournable dans le domaine. Nous avons choisi la première voie, c‟est à

dire de nous appuyer sur un code éléments finis orienté objet existant, FEMJava (voir

[Eyheramendy 03]), afin d‟y introduire les contributions liées aux nouvelles formulations

développées.

Dans ce contexte, les points clés de la génération automatique de code éléments finis sont

d‟une part, la réutilisation de code correspondant à la gestion globale des données dans le

code éléments finis, et d‟autre part, la particularisation du calcul des grandeurs élémentaires

développées sous forme symbolique. Ainsi, des extensions automatiques sous forme de

classes et de méthodes dans FEMJava permettent la particularisation du code sous la forme

d‟un précompilateur de formes symboliques. Le code créé sera pourvu du comportement

permettant de calculer les contributions élémentaires correspondant à la nouvelle formulation.

Nous cherchons ici à poser les premiers concepts d‟une approche homogène depuis les

équations du problème, jusqu‟aux simulations numériques. Le monde symbolique et le monde

numérique sont intégrés au sein d‟une même approche en Java.

Dans un premier temps, nous décrivons le code éléments finis orienté objet FEMJava. La

structuration de l‟application se prête au développement automatique de code de part sa

structuration en package thématiques. La présentation d‟une extension au code FEMJava

permet de poser les principes de programmation automatique. La mise en forme à objets est

présentée dans une seconde partie. Il est important de noter que ce schéma est complètement

général et peut être adapté à tous les langages de programmation et tout type codes éléments

finis, avec plus ou moins de facilité suivant la structure du code cible.

5.2 Un code éléments finis orienté objet en Java : FEMJava

FEMJava est une implémentation orientée objet de la méthode des éléments finis en Java

(voir par exemple [Eyheramendy 03, 04, 06, 07, 08]. Dans ce paragraphe, nous présentons les

principaux objets de FEMJava à travers de la description des packages de ce code éléments


75

finis. Nous ne rentrons pas dans les détails des classes, et nous nous contentons de décrire les

rôles des différents packages de manière extrêmement simplifiée.

5.2.1 Organisation des packages

FEMJava est organisé en packages thématiques (organisation hiérarchique). Le package

principal femjava regroupe l‟ensemble du code éléments finis. Il est montré par exemple dans

la Figure 53 Dans toutes les figures de ce paragraphe, les packages sont représentés en

italique et les classes en gras.

Les classes nécessaires à la description du maillage sont contenues dans le package

femjava.mesh (voir Figure 53).

femjava

mesh

Node

Vertex

Edge

ElementalField

ElementalNodalField

ElementalGaussPointFieldField

… assembling

AssembledOjbects

ElementalBlockMatrix

ElementalFullMAtrix delaunay

MeshDelaunay2DT1ViaHull

MeshDelaunayGenerator

MeshGeneratorL1 adaptivemesh

MetricFinitePointGenerator

Pointmeshenhacement elementalgeometry

Hexa1

Q0

Side …

meshing

MeshGenerator3D

MeshGeneratorQ1 …

Figure 53 Extrait du package femjava.mesh

Les classes Node, Vertex, Edge permettent de définir les nœuds, points, arêtes,… Les classes

ElementalNodalField et ElementalGaussPointField définissent les données locales des

champs scalaires, vectoriels ou tensoriels. Un mailleur 2D de Delaunay est programmé dans


76

les classes du sous-package femjava.mesh.delaunay. Dans le package

femjava.mesh.elementalgeometry, on trouvera les classes définissant les éléments

géométriques.

femjava

field

Field

ScalarField

VectorField

Vector2DField

Vector3DField

TensorField SubDomainField

SubDomainNodalField

SubDomainGaussPointField

operator

Gradient

Operator

Figure 54 Extrait du package fem.java.field

Dans le package femjava.field, on définit les champs (scalaire, vecteur ou tenseur). Les classes

contenues dans le package définissent les champs en tant que grandeurs globales (définies sur

le domaine). Cela permet de définir des champs aux nœuds ou des champs aux points de

gauss, suivant le type des points supports des données discrètes. Les données nodales ou aux

points de Gauss sont stockées dans des classes situées dans un autre package dévolu à cet

effet (voir Figure 54). Le sous-package femjava.field.operator contient les opérateurs

différentiels appliqués aux champs tel que le gradient et qui permet d‟effectuer des opérations

sur les champs (voir Figure 55).

femjava

field

Operator

Operator

Gradient

Figure 55 Extrait du package femjava.field.operator

Dans le package femjava.material, on définit les matériaux : élasticité linéaire, fluide

newtonien,…. Les classes définissant le comportement des matériaux sont contenues dans le

sous-package femjava.material.behaviors et les critères limites de lois de comportement sont

donnés dans femjava.material.behaviors.criterions (voir Figure 56).


77

femjava

material

LinearElastic

NewtionianFluid

MisesPerfectlyPlastic …

behaviors

Linear

linearElasticPlastic …

criterions

VonMises …

variable

InternalVariable integrators

RadialRetrnAlgorithm …

Figure 56 Extrait du package femjava.material

Le package femjava.formulation contient deux classes indispensables à la définition d‟une

nouvelle formulation. Ce sont les classes Element et Formulation (voir Figure 57). La

première permet de définir les données locales à un élément fini (nœuds,…), la seconde

contient des données plus générales ou globales (champs inconnues, champs annexes,

constructions des données élémentaires à partir des champs et caractéristiques de la

formulation). Dans la classe Element, on identifie l‟élément en déclarant toutes ses propriétés

telles que : la géométrie, la cinématique, le nombre d‟inconnues (nombre de champs), le

matériau de l‟élément, les charges volumiques….

femjava

formulation

Element

Formulation basic

LinearElasticity

NavierStokes

…

dim3

LinearElasticity3D

…

Figure 57 Le package femjava.formulation

La classe Element est une classe abstraite. Dans la classe Formulation, on initialise le

domaine et ses différents sous-domaines au travers de toutes les données de l‟élément et on


78

traite sa géométrie. C‟est également une classe abstraite, les formulations éléments finis sont

alors définies dans les sous-classes de celle-ci. Dans une sous-classe de la classe

Formulation. La classe Element sert de classe de base à une sous-classe définie au sein

d‟une sous-classe de la classe Formulation (mécanisme de classe interne). Ainsi, les données

locales (au niveau d‟un élément fini) et globales (les champs inconnus) nécessaires à la

définition d‟une formulation sont définies dans la sous classe de la classe Formulation ;

celle-ci contient une classe interne sous-classe de la classe Element qui elle contient toutes

les méthodes permettant de calculer les contributions élémentaires du problème éléments finis

Ces deux classes sont représentées dans la Figure 58.

Dans le package femjava.formulation.basic, on trouvera les classes représentant les

formulations tels que : élasticité linéaire, équations de Navier-Stokes incompressible… Les

formulations 3D sont isolées dans le sous-package femjava.formulation.basic.dim3.

Figure 58 Les classes Element et Formulation

Le package femjava.linearalgebra contient les classes d‟algèbre linéaire (voir Figure 59).

Dans ce package, la classe Matrix représente le concept de matrice. La classe FullMatrix est

l'abstraction d'une matrice pleine. Ces principales tâches sont liées au calcul algébrique

(addition, multiplication,…). La classe FullElementalMatrix hérite de sa superclasse les

tâches liées au calcul, et se particularise par la gestion de la numérotation liée à l'assemblage

de la contribution élémentaire pour le système linéaire global (position des inconnues dans le

système linéaire global).

Le sous-package femjava.linearalgebra.linearsolver contient les solveurs de systèmes

linéaires, solveur itératifs ou directs : méthode de Crout (stockage skyline), gradient conjugué,

Element

-geometry

-numberOfUnknownFields

-elementalField

-quadrature

-material

-kinematics

-bodyLoads

+setMaterial()

+getMaterial()

+setKinematics()

+…

Formulation

-material

-kinematics

+initialize(Domain)

+initialize(SubDomain)

+…


79

GMRES,… Le package femjava.linearalgebra.cuthillmckee contient des classes pour

renuméroter les inconnues du problème éléments finis.

femjava

linearalgebra

MathVector

Matrix

FullMatrix

… cuthillmckee

Edge

Vertex

…

linearsolver

IterativeSolver

ConjugateGradientMethod …

…

Figure 59 Extrait du package femjava.linearalgebra

Les formulations sont définies dans le package femjava.formulation.

5.2.2 Développement de nouveaux modèles éléments finis dans FEMJava

Dans ce paragraphe, nous allons décrire l‟implantation d‟un nouvel élément (nouveau modèle

d‟équation) dans FEMJava.

5.2.2.1 Principe d’intégration d’un nouvel élément dans FEMJava

La définition d'un nouveau modèle éléments finis nécessite la définition de deux types de

données:

Des données globales définies sur le domaine de calcul: définition des champs

inconnus et de champs annexes du problème

Des données locales élémentaires: calcul des matrices éléments finis en fonction de la

définition locale des champs, ...

Il est important de noter que la formulation entre dans un cadre forcément déjà implanté dans

le code élément finis, comme par exemple le fait d‟avoir un problème en statique ou

dynamique. La définition d'une formulation éléments finis se fait dans une sous-classe de la

classe Formulation. Au niveau de la classe, toutes les caractéristiques globales de la

formulation (champs inconnus, etc.…) sont données. Cette classe contient une classe interne

dans laquelle les éléments locaux sont définis : calcul des contributions élémentaires du

problème, matériau, schéma d‟intégration numérique,… Cette classe est une sous-classe de la


80

classe Element. A titre d‟exemple, nous donnons dans la suite un exemple d‟implémentation

de nouvelle formulation dans le code FEMJava. C‟est ce schéma que devra respecter notre

méthode de programmation automatique d‟un nouvel élément.

5.2.2.2 Application à un écoulement d’un fluide incompressible

On considère la formulation de Navier-Stokes pour un fluide newtonien incompressible. Les

équations du problème sont données dans la Figure 60. On cherche la vitesse et la pression

telles que les équations d‟équilibre, conditions de bords et condition d‟incompressibilité

soient respectées.

Figure 60 Forme forte des équations de Navier-Stokes pour un fluide incompressible

La formulation éléments finis utilisée pour développer le modèle numérique est une

formulation mixte (vitesse-pression), formulation de Galerkin stabilisée par ajout de termes de

type moindres-carrés (voir [Tezduyar 92]. Elle est donnée comme suit :

Considérons les espaces solution et test suivant :

Une formulation variationnelle mixte du problème est alors :

2

1

f

F

uu


81

où le paramètre de stabilisation est défini par :

L‟implémentation de la formulation est réalisée dans une classe appelée

NavierStokesFormulation. La classe est donnée dans la Figure 61. Elle appartient à un

package appelé femjava.fem.formulation.basics. La classe est une sous-classe de la classe

Formulation. La classe NavierStokesFormulation contient une classe interne statique

appelée NavierStokes. Les deux principales méthodes de la classe sont:

initialize() permet de décrire les champs inconnus, ici un champ scalaire, la pression,

et un champ vectoriel, la vitesse (2 inconnues la pression, champ scalaire, et la vitesse

champ vectoriel)

getElement() permet d'instancier au niveau local l'élément fini appelé NavierStokes

défini à partir d‟un certain nombre de données liées au problème : quadrature

numérique, matériau…

La classe NavierStokes est définie comme une classe interne statique (voir Annexe A Le

langage Java ou [Flanagan 06] pour plus de détails). Les instances de cette classe peuvent

ainsi avoir accès aux données de la formulation définie dans la classe

NavierStokesFormulation. C‟est une sous-classe de la classe Element qui permet de gérer

les données élémentaires: interpolations d'éléments finis,… La classe NavierStokes contient

toutes les méthodes permettant de calculer les contributions élémentaires liées à la

formulation éléments finis. Cette implémentation basée sur le concept de classe interne

permet au programmeur de définir les données au même niveau global et local pour le

problème, et ce au sein du même classe. Ceci garantit la cohérence entre élément et

formulation.


82

package femjava.fem.formulation.basics ;

public class LinearizedStabNSFormulation extends Formulation {

public String toString() { // … }

public Material defaultMaterial(){

return new NewtonianFluid(600.0, 1.0);}

public void initialize( Domain domain ) { // … }

public Element getElement( ElementalGeometry aGeom ,

Quadrature, aQuadrature , ElementalField[] flds,

int nb , Material m ) {

return new LinearizedStabilizedNavierStokes(aGeom,

aQuadrature, flds, nb, m, k); }

public static class LinearizedStabilizedNavierStokes

extends Element {

public LinearizedStabilizedNavierStokes(ElementalGeometry

geom, Quadrature aQuadraV, ElementalField[] flds,

int nb, Material m, Kinematics k) {

super(geom, aQuadraV, flds, nb, m, k);}

public void flow(TimeStep ts) {}

public FullMatrix computeT() {return null;}

public MathVector computeDelta() {return null;}

// … additional non-implemented abstract methods …}}

Figure 61 Implémentation d‟une classe interne pour définir une formulation

Pour conclure cette partie, on considère un écoulement d‟un fluide incompressible autour d‟un

cylindre (voir Figure 62). La vitesse est imposée en entrée d‟écoulement.

Figure 62 Le domaine avec le cylindre

Les données du matériau sont : la densité volumique , le module de

cisaillement . Le nombre de Reynolds pour cet écoulement est de 100 (rapporté au

diamètre). Dans la Figure 63, on montre les champs de vitesse et de pression à un instant

(régime non permanent, écoulement oscillant). Les résultats sont conformes à ceux que l‟on

trouve dans la littérature (période des oscillations).


83

Figure 63 Champ de vitesse et de température à .

5.3 Outil de génération automatique de code

Dans le paragraphe précédent, nous avons décrit l‟application éléments finis cible FEMJava

pour l‟intégration des données obtenues sous forme symbolique. Ce schéma d‟implémentation

est automatisé. Nous décrivons ici les principes de l‟implémentation automatique dans

FEMJava en nous appuyant sur l‟exemple de l‟élasticité linéaire. Nous commençons par

décrire l‟élaboration de la formulation sous forme symbolique sur l‟exemple de l‟élasticité

linéaire utilisé tout au long de ce document.

5.3.1 Développement d’une formulation variationnelle sous forme

symbolique

D‟un point de vue pratique, la formulation (instance de la classe Formulation, voir chapitre

précédent) est développée dans une classe, ici nommée de manière générique classe

FileGenerating, sous forme de script (voir Figure 64). Le compilateur Java sert

d‟interpréteur pour le langage de commande permettant de construire le problème

variationnel. Le langage de commande (langage de modélisation) s‟appuie sur les

développements de classes réalisés en Java. L‟objet formulation va ensuite automatiquement

générer un nouvel objet dans FEMJava qui décrit le problème et respecte la sémantique de

ce code cible.

FileGenerating

-formulation : Formulation

+main(String[] args)

Figure 64 Classe FileGenerating

Le développement de la formulation se fait dans une classe du type de celle présenté en

Figure 64. La classe FileGenerating conserve une instance de la classe Formulation qui

contient les éléments construits dans la fonction principale (fonction principale main(String[]


84

args)). Un extrait du contenu de la fonction main(String[] args) est présenté dans la Figure

65 et se décompose de la manière suivante :

On définit tous les champs de la formulation (code A).

On construit les termes de base de la formulation en appliquant les opérateurs différentiels

aux champs (code B).

On construit les tenseurs représentant l‟approximation éléments finis des champs (code C).

On effectue les opérations tensorielles entre champs (code D).

On crée une instance de la classe formulation (code E) qui contient en particulier tous les

termes de la formulation.

//------------------------------A----------------------------

1 Field N1 = new Field("N1", null) ; //




5 Field u = new Field("u", new int[] { 2 }, 2, true,

new Operator[] {N1, N2, N3, N4}); //

6 Field v = new Field("v", new int[] {2}, 2, true, null);//

7 Fraction lambda = new Fraction (173.0e6, 1.0);

8 Fraction mu = new Fraction (115.0e6, 1.0);

9 Tensor C = u.giveConstitutiveTensor(lambda,mu);//

10 Field f = new Field("f", new int[] {2}, 2, false,null);//

//-------------------------------B----------------------------

11 DifferentialOperator Eu =

new DifferentialOperator("E", u);//

12 DifferentialOperator Ev =

new DifferentialOperator("E", v);//

//-------------------------------C----------------------------

13 Tensor TEu = Eu.giveDiscretizationTensor(); //

14 Tensor TEv = Ev.giveDiscretizationTensor();//

15 Tensor Tv = v.giveDiscretizationTensor();//

//-------------------------------D----------------------------

16 Tensor TEuC = TEu.doubleContraction(C); //

17 Tensor TEuCTEv = TEuC.doubleContraction(TEv); //

18 Tensor fTv = moins.(f.simpleContraction(Tv)); //

//-------------------------------E----------------------------

19 Tensor[] leftSymbolicTensor= new Tensor[] {TEuCTEv};

20 Tensor[] rightSymbolicTensor= new Tensor[]{fTv};

21 Formulation g = new Formulation("LinearElasticity", new

Equation []{equation}, new Material (lambda, mu),

leftSymbolicTensor, rightSymbolicTensor);

Figure 65 Script décrivant une nouvelle formulation EF - Classe FileGenerating


85

Ce code est un script de définition des formulations. On peut considérer que cela représente

un langage de commande ou un langage de modélisation de formulation éléments finis.

Pour détailler les étapes de la construction de ce script, n considère le problème d‟élasticité

linéaire :

où le vecteur représente la solution du problème, le vecteur la fonction test et le vecteur

des charges. On note :

,

et où

et sont les coefficients de Lamé. Le développement de la formulation variationnelle de ce

problème est décrit en suivant les étapes de la Figure 65.

5.3.1.1 Définition des éléments de base du problème

La première étape du développement consiste à créer tous les champs de base et les constantes

du problème (Partie A, lignes 1 à .10 de Figure 65).

On définit tout d‟abord (ligne 1 à 4) les fonctions d‟interpolations nécessaires à

l‟interpolation du champ u sur l‟élément e de la manière suivante :

, où

les et sont les degrés de liberté correspondant au champ u sur un élément e et

sont les fonctions d‟interpolation. Les fonctions sont définies comme des instances de la

classe Field. Le nom du champ est donné sous forme d‟une chaine de caractères.

Les champs solution u et test v sont ensuite définis ligne 5 et 6. Le constructeur Field est

utilisé pour définir le champ. Les arguments du constructeur sont :

- le nom

- les dimensions de ce champ sont définies dans un tableau d‟entiers dont la dimension

est l‟ordre de tensorialité de ce champ, et dont le contenu donne les dimensions

- la nature du champ est un booléen qui précise si le champ est solution/test, ou non

- les fonctions d‟interpolation qui sont sous forme d‟un tableau d‟expression


86

Le champ est donné dans la Figure 66, c‟est un tenseur d‟ordre 1 (vecteur) et de dimension

2. Il est un champ solution donc il sera discrétisé sur les points d‟interpolations données. Les

fonctions d‟interpolation sont des polynômes de degré 1 dont on se contente de donner les

noms (fonctions d‟interpolation pour i=1 à 4).

Field

name: „„u‟‟ ;

dimension : new int []{2} ;

solutionOrTest: true ;

interpolationFunctions : new Operator[]{N1, N2 N3, N4} ;

Figure 66 Description du champ inconnue du problème d‟élasticité

Le champ est donné dans la Figure 67, c‟est un tenseur d‟ordre 1 (vecteur) et de dimension

2.

Field

name: „„v‟‟ ;


SolutionOrTest: true ;

Figure 67 Description du champ test (ou virtuel)

Le champ donné dans la Figure 68, est un vecteur de dimension 2. Il n‟est ni un champ

solution ni un champ test il représente le vecteur des charges.

Field

name: „„f‟‟ ;


SolutionOrTest: false ;

Figure 68 Le champ

Le champ est donné à partir de la méthode giveConstitutiveTensor(Fraction,Fraction) qui

passe en argument les coefficients de Lamé. Ces coefficients sont donnés dans les lignes 7 et

8 de la Figure 65.

Le tenseur constitutif est de la forme :


87

où est le symbole de Kronecker défini par :

Le symbole de Kronecker est donné par la méthode giveKronecker(int, int) de la classe Field

(voir Figure 69).

public Fraction giveKronecker(int n, int m) {

Fraction a = null;

if (n == m)

a = new Fraction(1.0, 1.0);

else

a = new Fraction(0.0, 1.0);

return a;

}

Figure 69 La méthode giveKronecker(int, int)

Des méthodes complémentaires sont implémentées pour définir le tenseur constitutif. La

méthode giveConstitutiveTensor(Fraction,Fraction) (voir donné Figure 70) consiste à

construire un tenseur d‟ordre 4. Cette méthode est programmée dans classe Field. Puis, on

définit ensuite chacun des éléments du tenseur grâce au symbole de Kronecker.

public Tensor giveConstitutiveTensor(Fraction lambda,Fraction mu)

{int n = this.getSpaceDimension();

Tensor C = new Tensor(new int[] { n, n, n, n });

Operator op1=null; Operator op2=null; Operator op3=null;

Operator op4=null; Operator op5=null; Operator op6=null;

for (int i = 1; i <= n; i++)

for (int j = 1; j <= n; j++)

for (int k = 1; k <= n; k++)

for (int l = 1; l <= n; l++) {

op1 = lambda.scalarProduct

(this.giveKronecker(i, j)).

scalarProduct

(this.giveKronecker(k, l));

op2 = this.giveKronecker(i, l).

scalarProduct

(this.giveKronecker(j, k));

op3 = this.giveKronecker(i, k).

scalarProduct

(this.giveKronecker(j, l));

op4 = op2.sum(op3);

op5 = op4.scalarProduct(mu);

op6 = op1.sum(op5);

C.setComponent(new int[]{i,j,k,l},op6);}

return C;}

Figure 70 La méthode giveConstitutiveTensor()


88

Les éléments de base du problème sont maintenant définis : le champ solution , le champ

test , le vecteur des charges et le tenseur constitutif .

5.3.1.2 Construction des termes de base du problème

Cette étape consiste à créer tous les termes du problème (Partie B, lignes 11 et 12 de la Figure

65). On crée une instance de la classe DifferentialOperator. Le nom de l‟opérateur appliqué

au champ et le champ lui même sont passés en argument du constructeur.

Le terme représente l‟opérateur différentiel tenseur de déformation appliqué au champ

vecteur (voir Figure 71).

DifferentialOperator

name: „„E‟‟ ;

field : u

Figure 71 L‟opérateur différentiel

Le tenseur obtenu est d‟ordre égale à l‟ordre de plus un, donc d‟ordre 2. La dimension de ce

tenseur est égale à la dimension du champ de base.

Afin de compléter la description, dans la Figure 42, on donne les composantes du tenseur

de dimension (construit à l‟aide de la méthode buildTensor()).

1,1=[du1_dx1];

1,2=[0.5*(du1_dx2 + du2_dx1)];

2,1=[0.5*(du2_dx1 + du1_dx2)];

2,2=[du2_dx2] ;

Figure 72 les composante de

De la même manière, on peut définir le terme avec ses composantes en fonction du

champ .

Tous les termes du problème variationnel sont décrits (opérateurs et champs).

5.3.1.3 Discrétisation des champs solutions et tests

Dans cette étape (Partie C, lignes 13 à 15 Figure 65), on construit les tenseurs représentant

l‟approximation par éléments finis des champs. Les termes concernés sont seulement ceux qui

dépendent du champ solution et du champ test. Le tenseur est construit à partir de la méthode

giveDiscretizationTensor(Field) qui passe en attribut le champ à discrétiser. On rappelle que

cette méthode est programmée dans la classe FEDiscretization.


89

Le tenseur représentant l‟approximation par éléments finis de terme est un tenseur

d‟ordre 3 (le tenseur est d‟ordre 2 et en ajoutant une unité) et de dimension .

Sa notation indicielle est où est l‟indice lié à la discrétisation. La dimension de cet

indice est (nombre de points d‟interpolation nombre de degrés de libertés).

Dans la Figure 73, on donne une partie des composantes du tenseur .

1,1,1=[dN1_dx1]; … ;1,1,8=[0.0];

1,2,1=[0.5*(dN1_dx2)]; … ;1,2,8=[0.5*(dN4_dx1)];

2,1,1=[0.5*(dN1_dx2)]; … ;2,1,8=[0.5*(dN4_dx1)];

2,2,1=[0.0]; … ;2,2,8=[dN4_dx2];

Figure 73 une partie des composantes du tenseur

De même pour le terme , on trouve un tenseur d‟ordre 3 et de dimension . Sa

notation indicielle est où est l‟indice lié à la discrétisation. La dimension de cet indice

est .

Pour le terme , , on trouve un tenseur d‟ordre 2 (le tenseur est d‟ordre 1) et de dimension

. Sa notation indicielle est .

Dans la Figure 74, on donne les composantes du tenseur .

1,1=[N1];1,2=[0];1,3=[N2];1,4=[0];

1,5=[N3];1,6=[0];1,7=[N4];1,8=[0];

2,1=[0];2,2=[N1];2,3=[0];2,4=[N2];

2,5=[0];2,6=[N3];2,7=[0];2,8=[N4];

Figure 74 Les composantes du tenseur .

5.3.1.4 Les opérations tensorielles pour construire les matrices élémentaires

Cette étape consiste à effectuer les opérations de contraction entre les tenseurs issus de la

discrétisation des termes du problème (Partie D, lignes 16 et 18 Figure 65).


est un tenseur d‟ordre 3 ( ) et de dimension

.




90

Le tenseur issu de la contraction simple entre le tenseur d‟ordre 1 et le tenseur d‟ordre 2


5.3.1.5 Définition de la formulation du problème

La formation est créée dans la partie E de la Figure 65. Une instance de la classe

Formulation est créée à partir :

- d‟un nom de formulation (chaine de caractères)

- des équations du problème

- matrices élémentaires calculées dans l‟étape précédente. Le premier représente le

membre droit de l‟équation et le deuxième le membre gauche.

- identification des caractéristiques du matériau choisi pour le problème est effectuée en

créant une instance de classe Material de FEMJava. Les attributs de cette classe sont

des nombres réels.

La formulation variationnelle du problème de l‟élasticité linéaire est donc définie. Les termes

de l‟équation sont : pour le membre de droite la contribution élémentaire du terme

, et pour le membre de gauche, la contribution élémentaire du terme

. Ces matrices élémentaires sont insérées dans des tableaux dans les lignes 19 et

20 de la Figure 65.

L‟instance de la classe Formulation est donnée dans la Figure 37.

Formulation

-name : „„LinearElasticityFormulation‟‟ ;

-equations : new Equation[]{equation} ;

-leftSymbolicTensor :new Tensor[]{ } ;

-rightSymbolicTensor : new Tensor[]{ } ;

-material : new Material(lambda, mu) ;

Figure 75 L‟instance de classe Formulation

Le script décrivant les tenseurs qui représentent la formulation variationnelle du problème est

donné dans la ligne 17 de la Figure 65.

5.3.2 Génération automatique de code à partir des formes symboliques

La formulation sait générer automatiquement le code source en Java qui s‟intègrera

automatiquement dans le code éléments finis FEMJava. La méthode

generateFEMJavaCode() dans la classe Formulation gère toutes les étapes de cette


91

intégration . Le résultat est une classe qui s‟intègre dans le code FEMJava. L‟exemple issu de

l‟élasticité linéaire est montré sous forme simplifiée dans la Figure 76 et la Figure 77. Un

extrait de la méthode generateFEMJavaCode() indiquant les princiaples actions pour générer

ce code est donné dans la Figure 78.

//----------------------------A-----------------------------

package femjava.fem.formulation.basics;

import femjava.algorithms.TimeStep;

import femjava.fem.formulation.Element;

…

import femjava.quadrature.GaussPoint;

import femjava.quadrature.Quadrature;

public class LinearElasticityFormulation extends Formulation

{…}

//-----------------------------B----------------------------

…

public void initialize(Domain domain) {

kinematics = new Kinematics.Dimension2();

Field[] fields = new Field[1];

fields[0] = domain.createAVector2DField(0);

domain.setFields(fields);

domain.setNumberOfUnknownFields(1);

this.finalizeDomain(domain);}

//-----------------------------C----------------------------

public Hashtable computeElementalMatrices(TimeStep ts){

Hashtable elementalMatrices = new Hashtable(3);

int matDimu = elementalField[0].computeNumberOfDofs();

ElementalFullMatrix Kwu = new

ElementalFullMatrix(matDimu,matDimu);

MathVector fu=new MathVector(matDimu);

ElementalBlockMatrix stiffness = new

ElementalBlockMatrix(Kwu,this);

ElementalBlockMatrix load = new

ElementalBlockMatrix(fu,this);

Kwu.equationNumbers(elementalField[0].

equationNumbers());

fu.equationNumbers(elementalField[0].

equationNumbers());

Figure 76 Code généré automatiquement et intégré au code FEMJava (première partie)


92

//-----------------------------D----------------------------

…

FullMatrix T1=compute1(shapeDer1, new double[]{lambda,mu});

FullMatrix T2=compute2(shapeDer1, new double[]{lambda,mu});

Kwu.plusSimpleProduct(T1,volume);

fu.plusProduct(T2,volume);

elementalMatrices.put("K",stiffness);

elementalMatrices.put("f",load);}

return elementalMatrices;}

//-----------------------------E----------------------------

private ElementalFullMatrix compute1(double[] shapeDer,

double[] constant){

ElementalFullMatrix answer =new ElementalFullMatrix(8,8);

double lambda = constant[0];

double mu = constant[1];

MathTensor C=this.giveTensorC(lambda, mu);

MathTensor E=this.giveTensorE(shapeDer);

MathTensor EC=E.doubleContract(C);

MathTensor T1=EC.doubleContract(E);

answer.atPut(1,1, T1.getComponent(new int []{1,1}));


…


return answer;}

Figure 77 Code généré automatiquement et intégré au code FEMJava (deuxième partie)

public String generateFEMJavaCode() {

String s = "";

s = this.importFEMJavaPackages()

+ this.giveInitializeMethod()

+ this.giveComputeElementalMatrices()

+ this.giveMatricesDimension()

+ this. giveElementalFullMatrices()

+ this. giveElementalBlockMatrices()

+ this. giveComputeShapeFunctionAt

+ this. giveFullMatrices()

+ this.giveAssemblyMatrices()

+ this.giveElementalMatrices()

return s;

}

Figure 78 La méthode generateFEMJavaCode()


93

5.3.2.1 Intégration de l’environnement symbolique dans l’environnement numérique

La classe LinearElasticityFormulation créée dans FEMJava, contient les objets de

l‟environnement numérique et symbolique. La première contrainte dans la méthode

generateFEMJavaCode() est d‟assurer la cohérence entre ces deux environnement.

Les packages de FEMJava présentés dans le paragraphe 5.2.1 permettent de définir tout type

de formulations. La méthode importFEMJavaPackages() importe les packages nécessaire

pour définir chaque nouvelle formulation et les classes nécessaires pour utiliser les méthodes

considérées.

Le package femjava.geometry permet de définir la géométrie du problème. Le matériau est

choisi dans le package femjava.materal. Les entités du maillage sont définies dans le package

femjava.mesh (voir Figure 79).

import femjava.geometry;

import femjava.material;

import femjava.mesh;

import femjava.linearalgebra;

import femjava.field;

…

Figure 79 Résultat de la méthode importFEMJavaPackages()

5.3.2.2 Particularisation les champs inconnus au problème

La méthode giveInitializeMethod() (voir Figure 80) intègre la méthode initialize(Domain)

dans la classe LinearElasticityFormulation. La méthode initialize(Domain), présentée dans

le paragraphe 5.2.2, associe les champs inconnus du problème de l‟environnement

symbolique à l‟environnement numérique. Elle crée les champs inconnus dans

l‟environnement numérique en se basant sur la formulation issu de l‟environnement

symbolique Cette méthode initialise le domaine considéré et définit les champs du problème.

Dans la Figure 80, on montre la méthode initialize(Domain) correspondant au problème de

l‟élasticité linéaire considéré (un seul champ inconnu, le déplacement vecteur de dimension

2).


94

public void initialize(Domain domain) {

kinematics = new Kinematics.Dimension2();

Field[] fields = new Field[1];

fields[0] = domain.createAVector2DField(0);

domain.setFields(fields);

domain.setNumberOfUnknownFields(1);

Figure 80 La méthode giveInitializeMethod()

5.3.2.3 Entêtes de méthodes pour le calcul des grandeurs élémentaires

La méthode giveComputeElementalMatrices() (voir Figure 81) intègre la méthode

computeElementalMatrices (TimeStep) dans la classe LinearElasticityFormulation. Dans la

Figure 81, on montre le détail de la méthode computeElementalMatrices (TimeStep)

correspondant au problème.

La méthode giveMatricesDimension() calcule les dimensions des matrices élémentaires du

problème (ligne 2). La méthode giveElementalFullMatrices() définit les formes des matrices

élémentaires (ligne 3 et 4). Cette méthode est détaillée dans la partie suivante. La méthode

giveElementalBlockMatrices() définit toutes les matrices blocs du problème : matrice de

rigidité, matrice de masse (pour les problèmes en dynamique), vecteur de charges,…(ligne 7

et 8) . Le calcul des grandeurs élémentaires s‟effectue à chaque point de Gauss (ligne 9, 10 et

11).

La méthode giveComputeShapeFunctionAt() calcule les fonctions d‟interpolations à un point

de Gauss. La méthode giveComputeShapeFunctionFirstOrderDerivativesAt() calcule les

dérivées premières des fonctions d‟interpolations. Le méthode giveFullMatrices() calcule les

matrice complètes du problème. Le détail de ce calcul est donné dans la partie suivante (ligne

15 et 16). La méthode giveAssemblyMatrices() assemble les matrices complètes dans les

blocks matrices convenable (ligne 17 et 18). La méthode giveElementalMatrices() renvoie les

matrices élémentaires (ligne 19 et 20).


95

public Hashtable computeElementalMatrices(TimeStep ts){

1 Hashtable elementalMatrices = new Hashtable(3);

2 int matDimu = elementalField[0].computeNumberOfDofs();

3 ElementalFullMatrix Kwu = new

ElementalFullMatrix(matDimu,matDimu);

4 MathVector fu=new MathVector(matDimu);

5 Kwu.equationNumbers(elementalField[0].equationNumbers());

6 fu.equationNumbers(elementalField[0].equationNumbers());

7 ElementalBlockMatrix stiffness = new

ElementalBlockMatrix(Kwu,this);

8 ElementalBlockMatrix load=new ElementalBlockMatrix(fu,this);

9 for (int i = 0; i < quadrature.numberOfGaussPoints(); i++) {

10 double[] gp = quadrature.point(i);

11 double w = quadrature.weight(i);

12 double [] shape = elementalField[0].getGeometry().

computeShapeFunctionsAt(gp);

13 double [] shapeDer = elementalField[0].getGeometry().

computeShapeFunctionsFirstOrderDerivativesAt(gp);

14 double volume = geometry.computeVolume(gp, w);

15 ElementalFullMatrix T1 = compute1(shapeDer,

new double[] {lambda,mu});

16 ElementalFullMatrix T2 = compute2(shape);

17 Kwu.plusSimpleProduct(T1,volume);

18 Kwu.plusSimpleProduct(T2,volume);}

19 elementalMatrices.put("K",stiffness);

20 elementalMatrices.put("f",load);}

return elementalMatrices;}

Figure 81 La méthode computeElementalMatricesMethod(TimeStep)

5.3.2.4 Calcul des contributions élémentaires aux points d’intégration

Dans cette partie, nous détaillons la génération du calcul des contributions élémentaires des

lignes 17 et 18 Figure 81. Nous considérons la contribution élémentaire du terme du membre

de gauche de l‟équation

. On a déjà donné la contribution élémentaire issu

de ce terme est

, où


96

Le code permettant l‟évaluation de H au point d‟intégration est donné Figure 82 (voir appel

de la méthode ligne 15 Figure 81).

private ElementalFullMatrix compute1(double[] shape,double[]

shapeDer,double[] constant){

1 ElementalFullMatrix answer=new ElementalFullMatrix(8,8);

2 double dN1_dx1=shapeDer[0];








10 double lambda = constant[1];

11 double mu = constant[2];

12 double T_1_1=dN1_dx1*(lambda + 2.0*mu)*dN1_dx1 +


(0.5*(dN1_dx2)* mu + 0.5*(dN1_dx2)* mu)*0.5*(dN1_dx2);

13 double T_8_8=(0.5*(dN4_dx1)*mu + 0.5*(dN4_dx1)*mu)*

0.5*(dN4_dx1)+(0.5*(dN4_dx1)*mu + 0.5*(dN4_dx1)*mu)*

0.5*(dN4_dx1) + dN4_dx2*( lambda + 2.0*mu)*dN4_dx2;

14 answer.atPut(1,2, T_1_1);

15 answer.atPut(1,2, T_8_8);

return answer;}

Figure 82 Extrait de la méthode compute1(double[],double[],double[]))

Cette contribution dépend des dérivées premières des fonctions d‟interpolation

et des

coefficients de Lamé et . Ces paramètres sont calculés à partir des variables intermédiaires

qui importent les valeurs de l‟environnement numérique. Les variables correspondant aux


97

dérivées d‟ordre premières des fonctions d‟interpolation sont données de la ligne 2 à la ligne 9

de la Figure 82. Ceux qui correspondent aux coefficients de Lamé sont donnés dans les lignes

10 et 11.

Les composantes de la matrice sont calculées dans lignes 12 et 13, et les résultats affectés

lignes 14 et 15.

5.3.3 Application numérique

Afin de valider la nouvelle formulation créée, on considère une plaque 2D en acier de

longueur , de la largeur . On applique des conditions de symétrie à sa

base et au côté gauche. On applique une charge répartie sur la partie supérieure verticale vers

le bas d‟intensité (voir Figure 83).

Les caractéristiques du matériau sont les suivantes : le coefficient de compressibilité

et le module de cisaillement .

Figure 83 Configuration géométrique du problème de plaque 2D

La plaque déformée est montrée dans la Figure 84. Les résultats sont comparés à ceux

trouvés à l‟aide de l‟application Abaqus (déplacements).


98

Figure 84 Déformation de la plaque

Chapitre 6 Approche symbolique pour formes variationnelles en grandes déformations

99

Chapitre 6 Approche symbolique pour

formes variationnelles en grandes

déformations

6.1 Approches symboliques et grandes transformations

Au cours des dernières années, d‟importants progrès ont été réalisés dans le domaine des

techniques formelles pour le calcul en grandes déformations. Les développements

mathématiques de ce type de formulation sont en général assez lourds. D‟un point logiciel, les

formulations permettant de traiter les grandes transformations sont en général assez difficiles

à mettre en œuvre suivant : matériau non linéaire, configurations de référence et

intermédiaires,… Cependant, on peut constater qu‟un certain nombre d‟éléments du

formalisme Lagrangien sont génériques et peuvent donc être systématisés. On peut noter dans

la littérature des travaux réalisés dans cette voie. Un des travaux les plus récents est celui de

[Korelc 10] qui présente une approche pour l‟automatisation des formulations éléments finis

des problèmes de contacts en grandes déformations. Une formulation de problèmes de

contacts est traitée de manière automatique sur une base de description symbolique de haut

niveau d'abstraction qui permet la dérivation du vecteur élémentaire résiduel et de la matrice


100

tangente. La procédure est implémentée dans un générateur de code automatique AceGen et

dans l‟environnement élément fini AceFEM [Korelc, 09].

Dans ce chapitre, nous allons brièvement rappeler les concepts permettant de traiter les

grandes transformations dans le contexte des principes du chapitre précédent. Le lecteur

pourra se rapporter à [Ibrahimbegovic 09], [Bathe 96] et [Liu 01] pour les développements

théoriques. Nous allons mettre en évidence les concepts de base afin de manipuler les formes

variationnelles en grandes transformations à partir d‟un exemple de formulation classique

issue de la littérature [Ibrahimbegovic 09].

6.2 Cinématique des grandes transformations

La cinématique a pour objectif essentiel d‟étudier le mouvement d‟un corps solide

déformable, c‟est-à-dire identifier les configurations successives occupées au cours de temps

. Dans le cadre choisi d‟un espace euclidien tridimensionnel, , la configuration d‟un corps

solide est considérée comme un ensemble de points matériels (voir Figure 85) où chaque

point matériel est identifié par son vecteur position dans un repère donné. A l‟instant , la

configuration initiale est définie par les coordonnées d‟un point matériel sur le corps, donné

par :

Figure 85 Configuration initiale et courante d‟un corps solide

Le mouvement d‟un corps solide peut alors être présenté comme l‟évolution de sa

configuration où la nouvelle position de chaque point matériel à l‟instant est définie par la

fonction , qui à associe la position du point matériel à l‟instant :


101

Dans une approche Lagrangienne, le mouvement d‟un corps solide déformable considéré

comme un ensemble de points matériels dont on connaît la position à l‟instant par

.

On définit le gradient de déformation , peut transformer un vecteur de la

configuration initiale en un autre vecteur dans la configuration courante.

6.3 Formulations forte et faible d’un problème en grandes

transformations

Dans le cas de la cinématique non linéaire, où les déplacements ne sont plus suffisamment

petits pour confondre la configuration initiale avec la configuration finale, les équations

d‟équilibre peuvent être exprimées dans différentes configurations.

6.3.1 Formulation forte

Les conditions d‟équilibre de forces et moments agissants sur un élément conduisent à :

Avec est le premier tenseur de contraintes de Piola-Kirchhoff, est la contrainte

de Cauchy, est les charges extérieures volumiques.

6.3.2 Formulation faible

La formulation faible des équations d‟équilibre peut être construite de la même manière qu‟en

petites perturbations (voir [Ibrahimbegovic 09] par exemple). La description matérielle du

principe des travaux virtuels en grands déplacements s‟écrit :

Avec :

est la déformation virtuelle de Green Lagrange, est le second

tenseur de contrainte de Piola-Kirchhoff, avec où est le tenseur

constitutif et

est le tenseur de Green-Lagrange, est le gradient de

déformation et et sont les charges extérieurs, respectivement volumiques et surfaciques.


102

6.3.3 Formulation faible linéarisée

On suppose que l‟on connaît déjà la configuration déformée à l‟instant , notée ,

on applique un petit déplacement incrémental sur un incrément qui conduit à

une nouvelle configuration déformée définie par . Par changement de coordonnées, il est

possible d‟introduire une représentation matérielle du déplacement incrémental :

Ceci permet de calculer facilement le gradient des déformations dans la configuration

déformée à partir de la description matérielle :

En raison du choix d‟un petit déplacement incrémental, superposé à la configuration

pour produire la configuration , on peut obtenir le même résultat en linéarisant le

gradient de déformation, qui peut être calculée en utilisant la dérivée directionnelle dans la

direction du déplacement incrémental :

Cette technique de linéarisation consistante peut être utilisée pour obtenir la représentation de

tout autre objet concernant la configuration . On peut calculer la formulation faible

linéarisée de la fonctionnelle :

On peut ainsi calculer :


103

En posant

, on peut donc écrire :

Donc le problème variationnel linéarisé devient :

avec :

et

.

On pose pour la suite:


104

6.4 Génération des contributions élémentaires

Dans ce paragraphe, nous appliquons le formalisme développé dans le Chapitre 3 pour

élaborer les contributions élémentaires. On considère la forme qui est une somme des

formes bilinéaires constituant la formulation. Le problème variationnel peut s‟écrire sous la

forme:

où :

avec :

On considère la première forme

. L‟approche adoptée au

Chapitre 3 conduit à écrire la forme discrétisée de par un tenseur d‟ordre 2 de la forme :

En identifiant les indices:

, , , ,


105

,

,

,

,

;

;

Donc:

De point de vue mathématique, est un vecteur de dimension 3 que l‟on souhaite interpoler

sur un hexaèdre (8 nœuds) : fonctions d‟interpolations . On note avec

les composantes du vecteur . On peut alors écrire l‟interpolation du champ sur

l‟élément e de la manière suivante :

, où les , et sont les degrés

de liberté au nœud correspondant au champ sur un élément . On pose :

, le

vecteur des inconnues nodales sur un élément de dimension 8. On peut alors écrire

l‟interpolation du vecteur sur l‟élément e tel que:


106

On considère maintenant le terme

. Son expression est une

matrice de dimension . On peut alors écrire la forme discrétisée correspondant à

où E est un tenseur d‟ordre 3 et de dimension . Sa notation

indicielle est où est l‟indice lié à la discrétisation. La dimension suivant cette direction

est (nombre de points d‟interpolation nombre de degrés de libertés).

Les composantes du tenseur permettant d‟écrire le produit matriciel sont

données dans la partie symbolique (voir le paragraphe 6.5).

De la même manière, on identifie la forme discrétisée correspondant à pour obtenir

. Puis en contractant tous les termes, on obtient le tenseur d‟ordre 2 de dimension

.

On considère maintenant la deuxième forme

. On peut écrire

pour sa forme discrétisée qui est un tenseur d‟ordre 2 de la forme :

On peut alors identifier les indices:

, , , ,

,

,

,

,

,


107

En développant l‟écriture de :

On obtient le tenseur de dimension :

où les crochets représentent les tenseurs provenant de la discrétisation de la forme.

L‟expression de la discrétisation des tenseurs pour ce terme suit le même principe que celui

du premier terme. Le calcul des tenseurs est représenté dans la partie symbolique.

Après avoir calculé les contributions élémentaires de tous les termes, nous pouvons calculer la

matrice élémentaire qui correspond à la forme en effectuant la somme de ces contributions:

On note que cette démarche s‟applique de la même manière à .

On montre ainsi que le formalisme développé dans le Chapitre 3 est suffisant pour décrire le

formalisme des formes variationnelles en grandes transformations. On se retrouve donc dans

la situation du chapitre précédent pour décrire les formes symboliques du problème. La même

infrastructure symbolique est donc utilisable telle quelle pour dériver ce type de problèmes.

Dans le paragraphe suivant, nous décrivons ce développement.

6.5 Description symbolique du problème

Nous avons montré dans le paragraphe précédent que de nouveaux objets n‟étaient pas

nécessaires pour développer des formulations en grandes transformations dans un cadre

variationnel. Dans ce paragraphe, nous décrivons la méthodologie pour décrire les termes qui

figurent dans la formulation (1) décrite dans le paragraphe 6.3.3. A partir de l‟approche à

objet décrite dans le Chapitre 4, nous exprimons les termes de cette formulation en fonction

des champs solution et test .


108

6.5.1 Définition des éléments de base de la formulation

On considère l‟application qui transforme la configuration initiale en configuration actuelle.

La fonction (position courante) peut être écrite comme la somme du vecteur position

initiale et du vecteur déplacement :

étant le champ déplacement et la position dans la configuration initiale. On définit le

gradient de déformation comme étant la dérivée de par rapport à :

Donc, le tenseur de déformation est exprimé en fonction du déplacement et de la matrice

d‟identité . La Figure 86 montre la construction de ce tenseur dans la méthode

giveDeformationTensor().

public Tensor giveDeformationTensor() {

DifferentialOperator gradu = new

DifferentialOperator("grad", this); //

Tensor I = this.giveIdendityTensor(); //

Tensor deformationTensor = gradu.sum(I); //

return deformationTensor; //

}

Figure 86 Construction du tenseur de déformation

A partir du tenseur de déformation , on peut définir le tenseur de Green-Lagrange

, avec est la transposée de et la matrice d‟identité.

public Tensor giveGreenLagrangeTensor () {

Tensor T1 = this.giveDeformationTensor();//

Tensor T2 = T1.transpose(); //

Tensor I = this.giveIdendityTensor();//

Tensor T1T2 = T1.simpleContraction(T2); //

Tensor answer = T1T2.minus(I); //

Fraction frac = new Fraction(1,2); //

Tensor GreenLagrange = frac.scalarProduct(answer);

return GreenLagrange; //

}

Figure 87 Construction du tenseur de Green-Lagrange


109

Dans la Figure 87, on définit les détails de la définition de ce tenseur (méthode

giveGreenLagrangeTensor()).

Le tenseur constitutif est un tenseur constant d‟ordre 4

qui dépend des coefficients du Lamé et .

Le tenseur de Piola-Kirchhoff 2 est obtenu en effectuant une double contraction entre le

tenseur constitutif et le tenseur de Green-Lagrange . La Figure 88 montre ces opérations

(méthode givePiolaKirchhoff2()).

public Tensor givePiolaKirchhoff2(Fraction lambda,

Fraction mu) {

Tensor T1 = this.giveConstitutiveTensor(lambda,mu);//

Tensor T2 = this.giveGreenLagrangeTensor();//

PiolaKirchhoff2= T1.doubleContraction(T2); //

return PiolaKirchhoff2;

}

Figure 88 Construction du tenseur (Piola-Kirchhoff 2)

6.5.2 Définitions des éléments discrétisés et contributions élémentaires

Dans cette partie, l‟outil de discrétisation automatique est appliqué sur les champs de la

formulation. Les contributions de chaque terme pour l‟élaboration des matrices élémentaires

sont générées afin d‟effectuer les opérations de contraction entre les tenseurs issus de la

discrétisation.

La définition des opérateurs différentiels , et en fonction des champs de base nous

permet de définir les termes et . La discrétisation du problème consiste à effectuer la

discrétisation de ces deux termes.

La construction de la forme issue de la discrétisation de la fonctionnelle

, est représentée dans la Figure 89 (méthode giveGammaTild()).


110

public Tensor giveGammaTild(Field w){



Tensor Tgradu =gradu.giveDiscretizationTensor();

Tensor TgraduT = Tgradu.transpose();//

DifferentialOperator gradw = new

DifferentialOperator("grad", w); //

Tensor Tgradw = gradw.

giveDiscretizationTensor();//

Tensor TgradwT = Tgradw.transpose();//

Tensor TgraduTgradwT = Tgradu.

simpleContraction(TgradwT);//

Tensor TgraduTTgradw = TgraduT.

simpleContraction(Tgradw);//

Tensor GammaTild1 =

TgradwTTgradu.sum(TgraduTTgradw);//

Fraction frac=new Fraction(1,2); //

Tensor GammaTild =

frac.scalar(GammaTild1); //

return GammaTild;//

}

Figure 89 Construction du tenseur de la discrétisation de

La discrétisation du tenseur conduit à un tenseur d‟ordre 4 et de dimension

, où la dimension 24 correspond au nombre de degrés de liberté issu de la discrétisation

du champs u. Dans le Tableau 1, on donne les composantes de ce tenseur sous forme

symbolique.

1,1,1,1=[0.5*(dN1_dx1*dN1_dx1+dN1_dx2*dN1_dx2 +dN1_dx3*dN1_dx3

+ dN1_dx1*dN1_dx1 + dN1_dx2*dN1_dx2 + dN1_dx3*dN1_dx3)];

1,1,3,24=[0.5*(dN1_dx1*dN8_dx1+dN1_dx2*dN8_dx2+ dN1_dx3*

dN8_dx3+ dN1_dx1*dN8_dx1+dN1_dx2*dN8_dx2 + dN1_dx3*dN8_dx3)];

3,1,1,1=[0.0];

3,24,3,24=[0.5*(dN8_dx1*dN8_dx1 + dN8_dx2*dN8_dx2 + dN8_dx3*

dN8_dx3+dN8_dx1*dN8_dx1 + dN8_dx2*dN8_dx2 + dN8_dx3*dN8_dx3)];

Tableau 1 Extrait de composantes du tenseur de la discrétisation de


111

De la même manière, la discrétisation de la fonctionnelle

est

donnée dans la Figure 90 (méthode giveGamma()).

public Tensor giveGamma(Tensor F){



Tensor Tgradu = gradu.

giveDiscretizationTensor();//

Tensor TgraduT = Tgradu.transpose();//

Tensor FT = F.transpose(); //

Tensor TgraduTF = TgraduT.

simpleContraction(F);//

Tensor TgraduFT = Tgradu.

simpleContraction(FT);//

Tensor Gamma1 = TgraduTF.sum(TgraduFT); //

Fraction frac = new Fraction(1,2);

Tensor Gamma = frac.scalar(Gamma1); //

return Gamma;//

}

Figure 90 Construction du tenseur de la discrétisation de

La discrétisation du tenseur donne un tenseur d‟ordre 3 et de dimension .

Dans le Tableau 2, on donne quelques composantes de ce tenseur.

1,1,1=[0.5*(dN1_dx1*(dN1_dx1*X_1 + dN2_dx1*X_4 + dN3_dx1*X_7 +

dN4_dx1*X_10 + dN5_dx1*X_13 + dN6_dx1*X_16 + dN7_dx1*X_19 + dN8_dx1*X_22 +

1.0)+ dN1_dx1*(dN1_dx1*X_1 + dN2_dx1*X_4 + dN3_dx1*X_7 + dN4_dx1*X_10 +

dN5_dx1*X_13 + dN6_dx1*X_16 + dN7_dx1*X_19 + dN8_dx1*X_22 + 1.0))];









Tableau 2 Extrait de composantes du tenseur de la discrétisation de


112

6.5.3 Génération de code

On a maintenant tous les éléments pour générer les contributions élémentaires pour les deux

formes élémentaires

et

Ces opérations sont implantées dans une classe annexe appelée LargeScaleFormulation

donnée dans la Figure 91. Dans cette classe, les méthodes de bases développées dans les deux

paragraphes précédents sont utilisées pour définir les termes de base de la formulation

variationnelle.

public class LargeTransformationFormulation {

public static void main(String[] args) throws IOException {

…

Field u = new Field("u", new int[] { 3 }, 3, true,

new Operator[] { new Field("N1", null),

new Field("N2", null), new Field("N3", null),

new Field("N4", null),new Field("N5", null),


new Field("N8", null)}); //

Field w = new Field("w", new int[] { 3 }, 3, true, null); //

Fraction lambda = new Fraction(1.0e6,1);

Fraction mu = new Fraction(0.3,1);

Tensor C = u.giveLagrangienElasticityTensor(lambda,mu); //

…

------------------------------A1-----------------------------

Tensor GuF = u.giveGamma(F);//

Tensor GwF = w.giveGamma(F);//

Tensor GuFC = GuF.doubleContraction(C);//

Tensor GuFCGwF = GuFC.doubleContraction(GwF);//

------------------------------A2-----------------------------

Tensor Gtilduw = u.giveGammaTild(w);//

Tensor GtilduwS = Gtilduw.doubleContraction(S);//

Tensor[] leftSymbolicTensor= new Tensor[] {GtilduwS,GuFCGwF };

Tensor[] rightSymbolicTensor= new Tensor[]{moinsGuFw,wb,wt};

-------------------------Formulation--------------------------

Formulation g = new Formulation("LargeScale",

new Equation[]{equation},new Material (lambda, mu),


Figure 91 Construction du problème

Les matrices élémentaires et

sont de dimension . A titre d‟exemple, on

montre la composante (1,1) du tenseur sous forme symbolique dans la Figure 92.


113

1,1=[0.5*(dN1_dx1*dN1_dx1 + dN1_dx2*dN1_dx2 + dN1_dx3*dN1_dx3 +

dN1_dx1*dN1_dx1 + dN1_dx2*dN1_dx2 + dN1_dx3*dN1_dx3)*((lambda +

2.0*mu)*0.5*((dN1_dx1*X_1 + dN2_dx1*X_4 + dN3_dx1*X_7 + dN4_dx1*X_10 +

dN5_dx1*X_13 + dN6_dx1*X_16 + dN7_dx1*X_19 + dN8_dx1*X_22 +

1.0)*(dN1_dx1*X_1 + dN2_dx1*X_4 + dN3_dx1*X_7 + dN4_dx1*X_10 +

dN5_dx1*X_13 + dN6_dx1*X_16 + dN7_dx1*X_19 + dN8_dx1*X_22 + 1.0) +

(dN1_dx2*X_1 + dN2_dx2*X_4 + dN3_dx2*X_7 + dN4_dx2*X_10 + dN5_dx2*X_13 +

dN6_dx2*X_16 + dN7_dx2*X_19 + dN8_dx2*X_22)*(dN1_dx2*X_1 + dN2_dx2*X_4 +


dN8_dx2*X_22) + (dN1_dx3*X_1 + dN2_dx3*X_4 + dN3_dx3*X_7 + dN4_dx3*X_10 +

dN5_dx3*X_13 + dN6_dx3*X_16 + dN7_dx3*X_19 + dN8_dx3*X_22)*(dN1_dx3*X_1 +


dN7_dx3*X_19 + dN8_dx3*X_22)-1.0) + lambda*0.5*((dN1_dx1*X_2 + dN2_dx1*X_5

+ dN3_dx1*X_8 + dN4_dx1*X_11 + dN5_dx1*X_14 + dN6_dx1*X_17 + dN7_dx1*X_20

+ dN8_dx1*X_23)*(dN1_dx1*X_2 + dN2_dx1*X_5 + dN3_dx1*X_8 + dN4_dx1*X_11 +

dN5_dx1*X_14 + dN6_dx1*X_17 + dN7_dx1*X_20 + dN8_dx1*X_23) + (dN1_dx2*X_2

+ dN2_dx2*X_5 + dN3_dx2*X_8 + dN4_dx2*X_11 + dN5_dx2*X_14 + dN6_dx2*X_17 +

dN7_dx2*X_20 + dN8_dx2*X_23 + 1.0)*(dN1_dx2*X_2 + dN2_dx2*X_5 +


dN8_dx2*X_23 + 1.0) + (dN1_dx3*X_2 + dN2_dx3*X_5 + dN3_dx3*X_8 +


dN8_dx3*X_23)*(dN1_dx3*X_2 + dN2_dx3*X_5 + dN3_dx3*X_8 + dN4_dx3*X_11 +

dN5_dx3*X_14 + dN6_dx3*X_17 + dN7_dx3*X_20 + dN8_dx3*X_23)-1.0) +

lambda*0.5*((dN1_dx1*X_3 + dN2_dx1*X_6 + dN3_dx1*X_9 + dN4_dx1*X_12 +

dN5_dx1*X_15 + dN6_dx1*X_18 + dN7_dx1*X_21 + dN8_dx1*X_24)*(dN1_dx1*X_3 +


dN7_dx1*X_21 + dN8_dx1*X_24) + (dN1_dx2*X_3 + dN2_dx2*X_6 + dN3_dx2*X_9 +


dN8_dx2*X_24)*(dN1_dx2*X_3 + dN2_dx2*X_6 + dN3_dx2*X_9 + dN4_dx2*X_12 +

dN5_dx2*X_15 + dN6_dx2*X_18 + dN7_dx2*X_21 + dN8_dx2*X_24) + (dN1_dx3*X_3

+ dN2_dx3*X_6 + dN3_dx3*X_9 + dN4_dx3*X_12 + dN5_dx3*X_15 + dN6_dx3*X_18 +

dN7_dx3*X_21 + dN8_dx3*X_24 + 1.0)*(dN1_dx3*X_3 + dN2_dx3*X_6 +


dN8_dx3*X_24 + 1.0)-1.0))];

Figure 92 Composante (1,1) du tenseur


114


On considère une poutre 3D en acier de longueur , de la largeur et de

hauteur . La poutre est encastrée à sa base. On applique une charge répartie sur sa

surface supérieure du haut vers le bas d‟intensité (voir Figure 93).

Figure 93 La géométrie du problème

Les caractéristiques du matériau sont les suivantes : le module de Young et

. La déformation de la poutre est représentée dans la Figure 94.

Figure 94 La configuration déformée

On vérifie que nous trouvons les mêmes déplacements qu‟avec le code Abaqus comme on le

montre dans les Figure 95 et Figure 96.


115

Figure 95 La géométrie du problème sur Abaqus

Figure 96 La configuration déformée sur Abaqus


116

Chapitre 7 Applications

117


7.1 Hyperélasticité

Noun nous proposons dans ce chapitre de décrire une formulation pour l‟hyperélasticité dans

le prolongement de celle développée dans le Chapitre 6.

7.1.1 Formulation Forte

Les matériaux hyperélastiques se caractérisent par l‟existence d‟un potentiel de déformation

ne dépendant que des invariants tensoriels , et , qui définit la loi de

comportement.

En adoptant le formalisme Lagrangien, les équations locales du problème s‟écrivent :

où est le second tenseur de Piola-Kirchhoff, est la densité, est le tenseur de

Cauchy-Green droit, est le tenseur gradient de transformation, est la transformation

de la configuration initiale à la configuration actuelle, est le vecteur normal à la surface


118

extérieure, est le déplacement. Le même formalisme qu‟en grandes transformations dans le

Chapitre 6 est adopté.

Dans ce modèle, le comportement du matériau ne dépend que du choix du potentiel de

déformation. Il existe un certain nombre de modèles hyperélastiques. Nous adoptons un

modèle Néo-Hookéen compressible. Le potentiel Néo-Hookéen est une fonction de deux

invariants tensoriels et et s‟écrit :

où est le module de cisaillement, est le module de compressibilité et .

7.1.2 Formulation faible

On peut exprimer le second tenseur de Piola-Kirchhoff comme suit :

D‟où :

On obtient enfin l‟expression du tenseur de Piola-Kirchhoff :

La formulation faible du problème s‟écrit :

Avec :


119

est la déformation virtuelle de Green Lagrange, est le

gradient de déformation et et sont les charges extérieurs, respectivement volumiques et

surfaciques.

7.1.3 Linéarisation de la forme faible

On applique la technique de linéarisation développée dans le Chapitre 6 et on obtient :

Le gradient des déformations dans la configuration déformée à partir de la description

matérielle :

où

est le tenseur d‟élasticité Lagrangien d‟ordre 4,

le tenseur de Piola-Kirchhoff 2 : ,

.

Le tenseur peut être décomposé comme suit :

On calcule :


120

Le terme se calcule alors comme suit :

Multiplions les deux membres de l‟équation par

D‟où finalement :

Donc la formulation faible linéarisé du problème devient :

7.1.4 Description symbolique du problème

7.1.4.1 Définition des éléments de base de la formulation

Dans ce paragraphe, on identifie tous les termes qui figurent dans la formulation. Les termes

,

et

sont définis dans le chapitre précédent. Il reste à définir et .


121

Pour calculer , on calcule

où est le tenseur de Cauchy-Green droit.

Le tenseur constitutif est alors:

La Figure 97 montre la construction de ce tenseur qui n‟est autre que la transposition

naturelle de la formule ci-dessus.

public Tensor giveLagrangienElasticityTensor(Fraction lambda,

Fraction mu){

Tensor C = this.giveCauchyGreenDroitTensor(); //

Tensor CmoinsUn = C.inverted(); //

Tensor CHyper1 = CmoinsUn.

tensorProduct(CmoinsUn); //

Tensor lambdaCHyper1 = lambda.

scalarProduct(CHyper1); //

Tensor F = this.giveDeformationTensor(); // Fraction j = new Fraction(F.giveDeterminant(),1);

double lnJ = j.giveLogarithm() ;

Operator cte = 2*mu.

minus(2*lambda.scalarProduct(lnJ));// Tensor J= CmoinsUn.tensorProductJ(CmoinsUn); //

Tensor CHyper2 = cte.scalarProduct(J); //

Tensor CHyper = lambdaCHyper1.

sum(CHyper2); //

return CHyper; } //

Figure 97 La construction du tenseur d‟élasticité Lagrangien

On calcule également le tenseur de Piola-Kirchhoff 2 (voir Figure 98), tenseur d‟ordre 2 et

qui s‟exprime par :

public Tensor givePiolaKirchhoff2HyperTensor(Fraction lambda,

Fraction mu) {

Tensor C = this.giveCauchyGreenDroitTensor(); //

Tensor CmoinsUn = C.inverted(); //

Tensor I = this.giveIdendityTensor(); //

Tensor S1 = I.minus(CmoinsUn); //

Tensor muS1 = mu.scalarProduct(S1); //

Tensor F = this.giveDeformationTensor(); //

Fraction j = new Fraction(F.giveDeterminant(),1);

double lnJ = j.giveLogarithm() ;

Operator lambdaLnJ = lambda.scalarProduct(lnJ); // Tensor S2 = lambdaLnJ.scalarProduct(CmoinsUn); //

Tensor PK2Hyper = muS1. sum(S2); //

return PK2Hyper; } //

Figure 98 La construction du tenseur de Piola-Kirchhoff 2


122

7.1.4.2 Définitions des éléments discrétisés et contributions élémentaires

Tous les termes sont calculés dans le paragraphe précédent. On applique l‟outil de

discrétisation automatique pour générer les contributions des termes contenant la solution et la

fonction test aux contributions élémentaires.

public class HyperElasticityFormulation {

public static void main(String[] args) throws IOException

{…

Field u = new Field("u", new int[] { 3 }, 3, true,

new Operator[] { new Field("N1", null),





Field w = new Field("w", new int[] { 3 }, 3, true, null); //

Fraction lambda = new Fraction (1.5e6,1); // Fraction mu = new Fraction (0.5e9,1); //

Tensor C=u giveLagrangienElasticityTensor(lambda, mu);//

Tensor F = u.giveDeformationTensor(); //

Tensor GuF = u.giveGamma(F);//

Tensor GwF = w.giveGamma(F);//

------------------------------A1-----------------------------

Tensor GuFC = GuF.doubleContraction(C);//

Tensor GuFCGwF = GuFC.doubleContraction(GwF);//

------------------------------A2-----------------------------

Tensor S = u. givePiolaKirchhoff2HyperTensor(lambda,mu); //

Tensor Gtilduw = u.giveGammaTild(w);//

Tensor GtilduwS = Gu.doubleContraction(S); //

Tensor[] leftSymbolicTensor= new Tensor[] {GuFCGwF,GtilduwS };

Tensor[] rightSymbolicTensor= new Tensor[]{moinsGuFw,wb,wt};

-------------------------Formulation--------------------------

Formulation g = new Formulation("HyperHelasticity", new

Equation[]{equation},new Material (lambda, mu),


Figure 99 Matrices élémentaires du problème d‟hyperélasticité

On calcule les matrices élémentaires qui correspondent à

et à

dans la classe

HyperElasticityFormulation (voir Figure 99).

Dans la Figure 99, le champ solution est défini et la définition de la discrétisation sur un

élément est donnée. De même pour le champ test.


123


On considère un bloc 3D en élastomère (considéré compressible) de dimension ,

et . Le bloc est fixé à sa base. On applique une charge répartie sur la

surface supérieure du gauche à droite d‟intensité (direction horizontale). Les

caractéristiques du matériau sont les suivantes : le coefficient de compressibilité

et le module de cisaillement .

La géométrie du problème ainsi que la configuration déformée est donnée dans la Figure 100.

Figure 100 Cisaillement d‟un bloc d‟élastomère dans FEMJava

On vérifie que l‟on trouve les mêmes déplacements de la partie supérieure du bloc que dans le

code Abaqus (voir les résultats dans la Figure 101).


124

Figure 101 Cisaillement d‟un bloc d‟élastomère dans Abaqus

7.2 Thermoélasticité

On considère un problème de mécanique en petite perturbation faisant intervenir la

température. Nous nous appuyons sur les deux principes de la thermodynamique pour dériver

les équations du problème.

7.2.1 Loi de comportement mécanique et loi de Fourier

La déformation totale peut se décomposer en partie mécanique et en partie thermique. On

suppose que :

où est la déformation mécanique et est la déformation thermique avec est le

coefficient de dilatation, est l‟évaluation de la température et est la matrice d‟identité.

L‟énergie libre est donnée sous la forme .

Le premier principe de la thermodynamique s‟écrit :

où est l‟énergie interne.

L‟énergie interne s‟écrit :


125

où est l‟entropie.

On remplace par son expression :

’ ù ’ x :

On suppose que est en fonction de et de ; . On peut écrire donc:

’ ù ’ x :

On peut exprimer le terme :

Donc on peut écrire :

Le second principe de la thermodynamique s‟écrit :

On peut développer :


126

’ ù

On introduit , la dissipation thermique de la forme :

Or

Donc on obtient :

Une condition suffisante afin de respecter l‟hypothèse , venant du second principe de la

thermodynamique, consiste à effectuer ce choix:

La première condition donne la loi de comportement mécanique :


127

La deuxième condition donne l‟expression de l‟entropie

La dernière condition donne la loi de Fourier :

Où est le coefficient de conductivité.

7.2.2 Equation de la chaleur

Suite aux hypothèses faites dans le paragraphe 7.2.1, l‟équation de la chaleur s‟écrit:

Or nous avons vu que l‟enpropie :

En dérivant par rapport au temps on trouve :

On remplace l‟expression dans l‟équation de la chaleur pour obtenir :

Soit ;

On pose la capacité calorifique :


128

’ ù ’ x :

On fait l‟hypothèse tel que se sépare en partie purement mécanique et en partie

purement thermique :

Enfin, on choisit une forme quadratique pour chacune des parties :

où est un paramètre constant et est défini par :

.

On peut donc en tirer les expressions des différentes dérivées partielles pour :

Cela permet d‟écrire l‟équation de la chaleur sous la forme finale:

Où la chaleur latente et la capacité calorifique définies telles que :


129

La relation du comportement mécanique s‟écrit :

7.2.3 Forme forte et forme faible

On résume les équations du problème de thermoélasticité comme suit :

La formulation faible issue de cette discrétisation est de la forme :

On discrétise en espace le domaine , et on approxime les champs de déplacements et et

les champs de température et respectivement par , , et .

Un schéma de résolution en temps de types différences finies peut être adopté pour le

problème.

7.2.4 Description symbolique du problème

7.2.4.1 Définition des éléments de base de la formulation

Les équations faibles du problème font intervenir les formes intégrales suivantes:


130

Le but est de générer les contributions élémentaires correspondantes à chacun de ces termes.

Nous avons déjà traité dans des exemples précédents la plupart de ces termes.

Il reste à créer :

Le tenseur

et le scalaire

.

Pour montrer l‟approche, le tenseur constitutif est ici choisi tel que les coefficients de Lamé

ne dépendent pas de la température :

Le calcul avec des expressions de dépendances en température des coefficients de Lamé ne

pose pas de problèmes supplémentaires. Les calculs seront simplement un peu plus lourds. De

même, on peut exprimer le terme qui représente la chaleur latente par et le terme

qui représente la capacité calorifique .

Dans ce contexte, la Figure 102 montre la construction du tenseur dans la méthode

giveLatentHeatTensor().

public Tensor giveChaleurLatenteTensor(Fraction alpha,

Fraction lambda, Fraction mu)

{

Tensor C = this.giveConstitutiveTensor(lambda,mu); //

Tensor I = this.giveIdendityTensor(); // Fraction moins = new Fraction(-1,1);

Fraction moinsalpha=moins.multiplication(alpha); // Tensor CI = C.doubleContraction(I);//

Tensor CC = moinsAlpha.scalarProduct(CI);//

return CC; //

}

Figure 102 Le construction du tenseur


131

7.2.4.2 Définitions des éléments discrétisés et contributions élémentaires associées

Dans les Figure 103 et Figure 104, on calcule les matrices élémentaires correspondant à la

discrétisation des termes ,

, ,

et . Le script

de la formulation est écrit dans la classe ThermoElasticityFormulation.

public class ThermoElasticityFormulation {

public static void main(String[] args) throws IOException

{…

Field u = new Field("u", new int[] { 3 }, 3, true, new

Operator[] { new Field("N1", null),





Field w = new Field("w", new int[] { 3 }, 3, true, null);//

Field T = new Field("T", null, 3, true, new

Operator[] { new Field("N1", null),





Field P = new Field("P", null, 3, true, null);//

Fraction alpha = new Fraction (10e-6,1); // Fraction rho = new Fraction (970,1); // Fraction c0 = new Fraction (2000,1); //

Fraction K = new Fraction (0.12,1); //K

Fraction lambda = new Fraction (1.510e6,1); // Fraction mu = new Fraction (0.510e9,1); //

Fraction moins = new Fraction (-1,1); // DifferentialOperator Eu=new DifferentialOperator("E",u);//

DifferentialOperator Ew=new DifferentialOperator("E",w);//

Tensor I=u.giveIdendityTensor();//

//----------------------------A1-----------------------------

Tensor TEu=Eu.giveDiscretizationTensor(); //

Tensor TEw=Ew.giveDiscretizationTensor(); //

Tensor C=u.giveConstitutiveTensor(lambda,mu); //

Tensor TEuC=TEu.doubleContraction(C);//

Tensor TEuCTEw=TEuC.doubleContraction(TEw);//

…

Figure 103 Classe ThermoElasticityFormulation (partie 1)


132

//----------------------------A2-----------------------------

Tensor TT=T.giveDiscretizationTensor(); //

Tensor TTI=TT.tensorProduct(I); //

Tensor moinsalphaTTI=

moins.scalarProduct(alpha.scalarProduct(TTI));//

Tensor alphaTTIC=alphaTTI.

doubleContraction(C);//

Tensor alphaTTICTEw= alphaTTIC.

doubleContraction(TEw);//

//----------------------------A3-----------------------------

Tensor TP = P.giveDiscretizationTensor(); //

Tensor TTTP = TT.tensorProduct(TP); //

Fraction Ce = moins.multiplication(c0); //

Tensor CeTTTP= Ce.scalarProduct(TTTP); //

Tensor rhoc0TTTP= rho.scalarProduct(c0TTTP); //

//----------------------------A4-----------------------------

DifferentialOperator gradT = new

DifferentialOperator("grad", T); //

DifferentialOperator gradP = new

DifferentialOperator("grad", P); //

Tensor TgradT=gradT.giveDiscretizationTensor(); //

Tensor KTgradT=K.scalarProduct(TgradT); //

Tensor TgradP=gradP.giveDiscretizationTensor(); //

Tensor KTgradTTgradP =KTgradT.

simpleContraction(TgradT);//

//----------------------------A5-----------------------------

Tensor CC=u.giveChaleurLatenteTensor(alpha,lambda, mu);//

Tensor CCTEu=CC.doubleContraction(TEu); //

Tensor rhoT=T.scalarProduct(rho); // Tensor rhoTCCTEu=rhoT.scalarProduct(CCTEu); //

Tensor rhoTCCTEuTP= rhoTCCTEu.tensorProduct(TP); //

…

Figure 104 Classe ThermoElasticityFormulation (partie 2)

7.2.5 Applications numériques

On considère un parallélépipède en élastomère de dimension , et

. Les données du matériau sont : la densité volumique , le

coefficient de dilatation , le coefficient de conductivité et la

capacité calorifique est .


133

Dans la Figure 105, on impose sur 3 faces du parallélépipède des conditions de symétrie en

déplacement. Pour la température, on impose sur deux faces opposées T=300 K.

Figure 105 Problème solution en température uniforme

Dans la Figure 106, la déformation du bloc est bien homogène et correspond au résultat

attendu en déplacement. La solution en température du problème est bien constante sur le

domaine T=300 K (la température est uniforme sur l figure, voir échelle). Ce calcul permet de

valider une partie du couplage thermique/mécanique dans la formulation.

Figure 106 Déformation et température pour le problème température uniforme

T=300 K

T=300 K


134

Dans la Figure 107, on garde les mêmes conditions de symétrie en déplacement et on impose

un gradient de température avec T=280 K sur la face fixé et T=300 K sur la face libre.

Figure 107 Problème avec solution en gradient de température

Dans la Figure 108, on montre la déformation du bloc et la solution en température. Les

déformations correspondent bien à celles engendrées par le gradient de température solution

du problème.

Figure 108 Déformée et solution du problème à gradient de température

Dans la Figure 109, on bloque les déplacements sur trois faces et on impose un gradient de

température avec T=280 K sur une face fixé et T=300 K sur la surface opposée.

T=280 K

T=300 K


135

Figure 109 Bloc encastré avec gradient de température

Dans la Figure 110, on montre la déformée du bloc et la solution en température. C‟est une

vérification qualitative de la formulation proposée.

Figure 110 Déformée et solution en température pour le bloc encastré

Enfin, nous vérifions qualitativement les solutions sur un problème instationnaire dans le

cadre du problème défini dans la Figure 109. Les conditions de bords sont T=290 K et T=300

K sur deux faces opposées et une température initiale de T=280 K est imposée sur reste du

domaine (voir Figure 111).

T=290 K

T=300 K

u=0

u=0

u=0


136

Figure 111 Evolution en déformée et température initiales et finales du problème

Chapitre 8 Conclusion

137


8.1 Sur une approche à objets généralisée pour la mécanique non

linéaire

8.1.1 Cadre général de l’approche

Dans ce travail, nous avons présenté une approche pour les dérivations des formulations

éléments finis dans un cadre variationnel au sein d‟un environnement symbolique. Cette

approche s‟apparente à une approche hybride symbolique/numérique. Nous avons développé

un environnement complet de calcul formel intégrant des outils d‟analyse tensorielle. Ce

travail s‟appuie sur un code éléments finis classique FEMJava (orienté objet en Java). De ce

point de vue, notre approche s‟apparente bien à celles proposées dans [Eyheramendy 95],

[Korelc 97] et [Logg 07]. La principale différence de fond avec ces approches précitées est le

choix fort d‟intégrer l‟environnement symbolique et environnement numérique au sein d‟un

seul et unique environnement, ici développé en Java. Concrètement, un nouveau package

contenant les outils de traitement de formes variationnelles est intégrée dans le package global

préexistant FEMJava. L‟objectif clair, visé par cela, est d‟évoluer au sein d‟un seul

environnement de programmation (un seul langage de programmation). Ainsi, toute la

problématique d‟un calcul éléments finis peut être appréhendée grâce à un seul concept de

programmation. Ce fait, d‟un point de vue gestion de code (maintenance, extension,

modification de code,…), est très un facteur extrêmement avantageux en terme d‟efficacité

que ce soit aux équipes de développement ou pour un programmeur devant intervenir dans


138

différentes parties de code. En calcul scientifique, lors de l‟élaboration d‟un nouveau solveur

pour un problème donné, on peut effectivement être amené à intégrer/modifier n‟importe

lequel des algorithmes d‟un code. Ce point est différent des approches [Eyheramendy 95],

[Korelc 97] et [Logg 07] pour lesquels plusieurs application et langages de programmation (à

différents niveaux) sont utilisés.

8.1.2 Une approche tensorielle générique pour la dérivation de modèles

éléments finis

Dans ce travail, l‟approche symbolique de dérivation de modèles éléments finis est basée sur

l‟analyse tensorielle. Le formalisme utilisé est un formalisme mathématique assez classique.

Celui-ci est exploité de façon à intégrer ce formalisme dans la méthode des éléments finis. Le

champ auquel est appliqué un opérateur différentiel est écrit sous forme discrète sur un

élément. Ici, seule la définition mathématique de l‟interpolation d‟un champ est utilisée come

base au développement des formes discrètes. Cela permet de s‟affranchir des écritures

classiques des modèles éléments finis basées sur la notation Voigt, et s‟appuyer simplement

sur les opérations tensorielles classiques entre champs (produit tensoriel, produit simplement

contracté, produit doublement contacté,…). Le caractère très générique de l‟approche pour

obtenir les formes discrétisées de termes variationnels quelconques est, à notre connaissance,

original dans le contexte donné.

A titre de comparaison, dans [Eyheramendy 97], une simple interprétation des termes de la

forme variationnelle permet de construire les formes élémentaires. On retrouve par exemple

dans la description des contributions élémentaires la notation classique de Voigt conduisant à

la forme connue

(voir par exemple [Hughes 87]), tout comme dans [Korelc 97].

Ces deux approches sont très générales mais nécessitent une étude de cas sur les termes

variationnels traités. Le formalisme proposé par [Logg 07] s‟appuie sur la même base de

description tensorielle que celle que nous avons développée pour générer les contributions

élémentaires. Par contre, il ne s‟applique que dans le cas particulier où la fonction de

changement de repère local-global est une fonction affine. Le tenseur élémentaire est

décomposé en deux parties distinctes, la première dépendant des fonctions d‟interpolation et

la deuxième dépendant des composantes des deux repères. Par exemple, pour l‟utilisation

d‟interpolation de type Hermite, le formalisme de [Logg 07] doit être adapté. Cela nécessite

donc une étude de cas, contrairement à notre approche.


139

8.1.3 Un environnement à objets de calcul formel pour la dérivation de

formulations éléments finis

Un environnement à objets a été créé pour gérer d‟une part la représentation des formes

variationnelles et des formes discrètes (matrices élémentaire), et d‟autre part, tous les calculs

sous forme symbolique. La représentation des expressions est basée sur une structure

arborescente s‟articulant autour d‟opérateurs unaires et binaires. Le développement clé de

l‟environnement consiste à avoir crée un objet tenseur qui, d‟une part, permet de décrire les

éléments de réduction d‟un tenseur et, d‟autre part implémente les différentes méthodes de

contraction d‟indices et les opérations entre les tenseurs. Il est important de préciser que ce

sont les opérations qui permettent d‟obtenir les formes élémentaires, en tenant compte des

éléments liés à la discrétisation des champs des termes de la forme variationnelle (classe

FEDiscretization). L‟implémentation proposée est basée sur les éléments d‟analyse

tensorielle énoncés dans ce travail. On hérite donc des aspects génériques de ceux-ci.

L‟implémentation est réalisée en Java dans le but d‟être intégrée naturellement dans un code

éléments existant FEMJava lui aussi développé en Java.

Cette approche est assez similaire à celle développée dans [Eyheramendy 95, 97] mais offre

un caractère beaucoup plus général, d‟une part par la structure de représentation des

expressions beaucoup plus générale, et d‟autre part, par le caractère générique du formalisme

mathématique basé sur l‟analyse tensorielle. Dans [Korelc 97], l‟outil Mathematica est utilisé

pour réaliser les manipulations, ce qui donne à l‟environnement une puissance importante en

termes de manipulations symbolique. Il est structuré par des fonctions (et non des objets)

permettant d‟obtenir les formes élémentaires. De plus, il s‟appuie sur des environnements

basés sur plusieurs langages de programmation qui rend plus difficile la prise en main de

l‟environnement global. Dans [Logg 07], l‟implémentation pour générer les formes

symboliques est réalisée en C++, mais se base sur une décomposition différente de

l‟élaboration des matrices élémentaires qui ne permettent pas forcément le traitement de

toutes les formes d‟interpolation (provient du choix de faire intervenir les espaces

d‟interpolation dans la forme matricielle).

8.1.4 Une approche de génération automatique de code éléments finis

Enfin, les contributions élémentaires obtenues symboliquement dans l‟environnement de

calcul formel sont intégrées dans un code éléments finis classique existant, FEMJava. Nous

avons fait le choix de conserver une approche logicielle classique (code de calcul) pour

effectuer les calculs numériques. Ce type d‟approche a aujourd‟hui fait preuve de son


140

efficacité pour la simulation de problèmes. On y adjoint donc un environnement de calcul

formel. L‟environnement symbolique pour le traitement des formulations éléments finis et

l‟outil de codage automatique sont implantés en Java. Le code cible de simulation numérique

est également implanté en Java. Le lien entre les deux mondes symboliques et numériques, est

obtenue grâce à un concept orienté objet pour la programmation automatique des formes

élémentaires. Le principe de programmation automatique est simple puisqu‟il suffit de

générer les entêtes du fichier correspondant à la nouvelle formulation et d‟intégrer les

expressions correspondant aux matrices élémentaires. Il faut noter que c‟est à ce niveau que

sont choisis les quadratures et les fonctions de changement de repères local/global. Cela n‟est

pas un inconvénient en soit, car il est nécessaire que le code cible doit pouvoir gérer les types

de géométrie et caractéristiques associées. Une généralisation de ces deux aspects sous forme

symbolique pourrait être facilement envisagée, par exemple intégrée au niveau de l‟objet

Integral (là ou se gère le schéma d‟intégration numérique et le changement de variable

local/global de la méthode des éléments finis).

La principale différence avec les approches similaires ([Eyheramendy 95], [Korelc 97] et

[Logg 07]) est le fait d‟avoir intégré naturellement les modèles équationnels au sein d‟un

environnement numérique, les deux étant basés sur un même paradigme de programmation.

Cela permet d‟envisager à plus long terme une prise en charge de composant du code de

calcul numérique directement dans la partie symbolique, et d‟y conserver la cohérence globale

de la description mathématique du problème éléments finis (éléments géométriques, schémas

de quadrature,…).

8.1.5 Sur la validation de l’approche proposée

Au cours de ce travail, l‟approche a été validée sur de nombreux problèmes classiques de la

mécanique, problèmes linéaires et non linéaires (élasticité, équation de la chaleur,...), et dans

le cadre de problèmes plus originaux, traitement de formulations variationnelles mixtes

(thermomécanique, Navier-Stokes,…) et grandes transformations (élasticité en grandes

transformations, hyperélasticité,…). Les applications choisies ici tentent de couvrir les

grandes tendances actuelles des problèmes numériques en mécanique des milieux continus :

grandes transformations, problèmes multiphysiques. On peut noter que le cadre Lagrangien a

pu être abordé grâce aux opérateurs différentiels classiques qui par combinaison ont permis de

construire les différents termes variationnels de la formulation Lagrangienne totale utilisée ici.

A la vue de la diversité des problèmes de la mécanique traités ici, on peut affirmer que


141

l‟application du paradigme objet d‟une part, et le choix du formalisme tensoriel d‟autre part,

sont pertinents pour couvrir une large gamme de problèmes. Bien qu‟il soit évidemment

impossible de développer un environnement capable de résoudre tout type de problème, nous

avons montré qu‟il était possible de construire un environnement structuré dans lequel il est

facile d‟intégrer des concepts nouveaux (nouveaux opérateurs différentiels, nouvelles

méthodes de discrétisation,…). Tous les points d‟entrées existent dans l‟approche proposée.

Contrairement aux approches proposées dans la même problématique ([Eyheramendy 97],

[Korelc 97] et [Logg 07], il semble évident que le traitement d‟un nouveau problème

n‟entraine que l‟ajout de nouvelle fonctionnalités aux objets existants, et ajout de nouvelles

classes, et cela dans la même cohérence globale. Le fait qui permette d‟affirmer cela est qu‟à

partir des concepts énoncés sur le cas simple de l‟élasticité linéaire, nous avons pu développer

des modèles éléments finis pour des problèmes plus complexes en enrichissant

l‟environnement symbolique des opérateurs adéquat (additions, produits ou opérateur

différentiel) formulation mixte pour les équation de Navier-Stokes en incompressible (voir

Annexe 2), formulation lagrangienne totale avec un modèle matériau hyperélastique,

problème de thermomécanique couplé. De même, le schéma de programmation automatique

dépend du code cible mais les concepts liés à la programmation automatique de code éléments

finis sont généraux et pourraient être enrichi de technique de gestions de variable temporaires

plus élaborés que l‟on peut trouver dans [Korelc 97] ou [Logg 07], voire de techniques de

différentiation automatique comme proposé dans [Korelc 97].

8.2 Pour aller plus loin dans ce travail

L‟ensemble de la démarche est intégré au sein d‟une approche unique basée sur le paradigme

orienté objet. Un formalisme mathématique a été développé afin de gérer automatiquement

l‟élaboration des matrices élémentaires d‟une formulation éléments finis pour un problème

physique donné. Ce formalisme est basé sur du calcul d‟algèbre tensorielle étendu à la

description de schémas de discrétisation éléments finis. Au regard de la littérature du domaine

(principalement [Eyheramendy 97], [Korelc 97] et [Logg 07], et au vue de l‟expérience, il

sera vraisemblablement nécessaire d‟enrichir à court terme l‟environnement de soit de

techniques de simplification d‟expression (notion que l‟environnement Mathematica sur

lequel s‟appuie [Korelc 97] intègre naturellement ou qui a été intégré dans [Logg 07]), soit de

techniques de différentiation automatique afin de simplifier le codes généré automatiquement.

Un système de gestion automatique de variables intermédiaires lors de la génération de code


142

plus performant que celui développé à ce jour sera aussi nécessaire. Il faut cependant noter

que pour les problèmes traités cela n‟a pas posé de problèmes de taille de code généré

(croissance rapide des expressions due à l‟approche tensorielle adoptée).

Dans notre approche, le développement du modèle éléments finis est basé sur l‟écriture d‟un

script en Java dans lequel sont décrits d‟une part la forme continue du problème, et d‟autre

part, les formes discrètes de chacun des termes variationnel. Il sera indispensable à court

terme de proposer une automatisation de la construction des formes élémentaires à partir de la

forme continue. La conséquence de cela sera une simplification du langage de modélisation,

dans lequel il ne sera plus nécessaire de décrire les produits de termes de la forme élémentaire

mais de désigner les termes à discrétiser dans un produit de termes. C‟est vraissemblablement

l‟élément à predre en compte le plus rapidement pour un développement futur.

Dans le même cadre de problématique, les problèmes non linéaires sont linéarisés dans une

phase d‟étude préliminaire (sur papier) avant leur introduction dans le script du problème. Il a

été montré dans le passé que des techniques de linéarisation basées sur la dérivée

directionnelle (voir [Eyheramendy 97] ou [Logg 07]) pouvaient être très naturellement

intégrées dans un environnement de calcul formel.

L‟approche que nous avons présentée ici traite la discrétisation par éléments finis des

équations de la physique dans un cadre variationnel. L‟extension à d‟autres types de méthodes

de discrétisation semble naturelle dans l‟approche mais reste à démontrer, par exemple par

adjonction de fonctions non polynomiales sur des supports variables (méthodes spectrales,

méthodes isogéométriques,…).

Par contre, un aspect essentiel des outils de simulation reste encore le traitement

algorithmique généré pour la résolution d‟un problème donné. L‟extension de cette réflexion à

la description algorithmique d‟un code de calcul devrait permettre l‟intégration au niveau

symbolique de nouveaux schémas, comme par exemple des schémas d'intégration en temps,

stratégies de calcul parallèle,… Ce sont là des ingrédients importants pour le calcul en haute

performance, lui-même incontournable aujourd‟hui d‟un point de vue industriel ou en

recherche. On peut envisager dans ce sens deux voies de recherche complémentaires, la

première basée sur une approche mathématique de la description des algorithmes, la seconde

sur une approche plus informatique de classification d‟algorithmes.


143

Enfin, l‟intégration de la dérivation de formes élémentaire à partir des équations d‟un

problème physique dans un environnement de type GUI (Graphical User Interface) comme

dans [Eyheramendy 97] permettrait de gagner en convivialité pour l‟élaboration de code de

calcul.

8.3 Les outils de calcul de demain

Dans le contexte global et intégré présenté dans ce travail, on ne s‟intéresse plus seulement à

créer une application particulière pour un problème donné, mais à concevoir de manière

globale un environnement intégrant des méthodes de traitement de problèmes issus de la

physique et des classes de méthodes de résolution algorithmique. Equations d‟un problème

variationnel, dérivations formelles des contributions élémentaires et simulations numériques

se retrouvent donc réunis au sein d‟un seul et unique environnement. Nous pensons que cela

préfigure la forme des outils de simulation que l‟on trouvera dans les années à venir. Ce type

d‟outils ou d‟approches doit tendre à libérer le mécanicien ou le numéricien des contingences

liés à l‟outil logiciel pour lui permettre de se concentrer sur son cœur de métier : la simulation

d‟un problème physique donné.


144

Annexe 1

145

Annexe 1 Java, un langage à objets

pour le calcul scientifique

Dans cette partie, on présente les principaux concepts du langage orienté objet Java. On met

en évidence les principales fonctionnalités de ce langage dans le contexte du calcul

scientifique.

1.1 Quelques généralités sur le langage Java

Java a été conçu pour être un langage sécurisé. Les caractéristiques de sécurité qui attirent

l‟attention sur Java concernent la possibilité d‟avoir un programme dynamiquement portable.

Il fournit plusieurs couches de protections contre les codes dangereux tels que les virus. Ces

mécanismes sont les fondations d‟une politique de sécurité de haut niveau qui autorise ou non

certaines activités en fonctions des tâches qu‟une application doit effectuer.

1.1.1 Sécurité de typage

Dans un langage statiquement typé comme C ou C++, les types de données sont

physiquement stockés dans le binaire lors de la compilation du code source. Le compilateur

peut identifier un bon nombre d‟erreurs avant exécution du code. Mais, ces langages ne

permettent pas la gestion de constructions de haut niveau. Cela rend impossible qu‟une

application intègre de nouveaux types de données en toute sécurité lors de l‟exécution.

Annexe 1

146

Un langage dynamiquement typé comme Smalltalk possède un système d‟exécution qui gère

le type des objets et effectue la vérification de type lors de l‟exécution de l‟application. Ce

type de langage permet des comportements plus complexes mais il est plus lent.

Java est un langage statiquement typé et à liaisons dynamiques. Chaque objet possède un type

bien défini qui est connu lors de la compilation. Le système d‟exécution de Java garde une

trace de tous les objets, et permet de connaître leurs types et leurs relations durant l‟exécution.

Les changements d‟un type d‟objet à un autre sont vérifiés par le système. Il est possible

d‟utiliser de nouveaux types d‟objet chargés dynamiquement.

1.1.2 Gestionnaire de Sécurité

Java possède un gestionnaire de sécurité, c‟est un objet qui peut être installé par une

application pour restreindre l‟accès aux ressources du système. Il est surtout utile pour des

applications qui lancent du code douteux lors de leur fonctionnement courant. L‟intégrité d‟un

gestionnaire de sécurité repose sur la protection offerte par les couches de bas niveau du

modèle de sécurité de Java. Cela peut être intéressant dans le cadre de calcul distribué.

1.1.3 Portabilité

Dans la plupart des langages, on dit qu‟un programme est portable si un même code source

peut être exploité dans des environnements différents moyennant une nouvelle compilation.

Dans ce langage, la portabilité va plus loin. En effet, la compilation d‟un code source produit

non pas des instructions machine, mais un code intermédiaire formé de byte code. D‟une part,

ce code est exactement le même quelque soit le compilateur et l‟environnement concernés.

D‟autre part, ce byte code est exécutable dans toute implémentation disposant du logiciel

d‟interprétation nommé machine virtuelle.

1.1.4 Machine virtuelle

Java est à la fois un langage compilé et interprété. Le code source est transformé en de

simples instructions binaires, qui ressemblent à du code machine. Le code est compilé en un

code universel, c'est-à-dire en instructions destinées à une machine virtuelle.

L‟application est exécutée par un interpréteur de Java. Cet outil fonctionne comme un

environnement virtuel sécurisé (Figure 112). Il peut être implémenté en n‟importe quel

langage.

Annexe 1

147

Figure 112 Exemple du fonctionnement de JVM

Java a été conçu pour que l‟implémentation du programme d‟exécution puisse optimiser

l‟exécution en compilant le pseudo-code en code machine. Cela s‟appelle une compilation

juste à temps (JIT). Le problème avec cette compilation est que l‟optimisation de code prend

du temps ce qui est un handicap pour obtenir de bonnes performances. Une nouvelle machine

virtuelle nommé Hotspot résout le problème en exploitant ce qu‟on appelle la compilation

adaptive. Ce morceau de code peut ne représenter qu‟un faible pourcentage du programma

total, mais son fonctionnement est déterminant pour sa performance globale.

1.1.5 Gestion dynamique de mémoire

La plupart des grandes différences entre Java et C ou C++ réside dans la gestion de mémoire

qui est complètement automatisée. Java intègre un système de ramasse-miettes basé sur une

gestion efficace des références. Ces caractéristiques éliminent les problèmes qui mettent en

danger la sécurité, la portabilité et l‟optimisation de code.

Les tableaux en Java sont des objets de première classe. Ils peuvent être alloués

dynamiquement et affectés comme tout autre objet. Les tableaux connaissent leur propre taille

Annexe 1

148

et leur type. Ils possèdent des relations d‟héritage bien définies et fondées sur les relations de

leur type de bases.

1.2 Java : un langage orientée objet pour la simulation en

mécanique

Java est un langage simple qui ne possède aucune ambiguïté de syntaxe. Il possède les

concepts classiques de la programmation orienté objet auxquels s‟ajoutent des structures de

langages propres.

1.2.1 Objets et Classes

La programmation orienté objet contribue à la fiabilité des logiciels et facilite la réutilisation

de code existant. Elle facilite en outre le développement à grande échelle d‟applications de

grandes tailles (plusieurs centaine de milliers de lignes de code). En programmation orientée

objet, chaque objet associe des données et des méthodes agissant exclusivement sur les

données de l‟objet. Il n‟est pas possible d‟agir directement sur les données d‟un objet ; il est

nécessaire de passer par ses méthodes. Ce concept est nommé encapsulation des données. Il

présente un intérêt manifeste en matière de qualité de logiciel. L‟encapsulation des données

facilite la maintenance : une modification éventuelle de la structure des données d‟un objet

n‟a d‟incidence que sur l‟objet lui-même. En Java, l‟encapsulation se trouve naturellement

induite par la syntaxe du langage mais elle n‟est pas absolue (comme dans bien des langages

orientés objet). L‟unité fondamentale de code Java est la classe qui est la description d‟un

ensemble d‟objets ayant une structure de donnée commune et disposant des mêmes méthodes.

A titre d‟illustration, on considère la classe Node présentée dans la Figure 113. Les champs

de la classe sont divisés en attributs et méthodes. Les attributs sont les coordonnées, les

valeurs de degré de liberté et les numéros équations dans un système linéaire donné pour les

attributs. Les méthodes nous permettent de gérer les coordonnées et différents types de

collecte. Sur la même figure, on donne deux exemples d'instances de la classe Node, appelé

ici node 23 et node 345. Chacun d'entre eux a ses propres attributs, c'est-à-dire ici, pour le

nœud 23, les coordonnées sont {2.5; 1.8}, les valeurs de degré de liberté (ici un champ de

vecteur à 2 dimensions) sont {0.025; 1.25}, et les positions dans le système linéaire (positions

des inconnues dans un vecteur) {732; 478}.

Annexe 1

149

Node 23

coordinates : [ 2.5; 1.8]

dofsValues : [0.025; 1.25]

equationNumber:[732;478]

Node 345

coordinates : [14.5; 42.3]

dofsValues : [2.567; 2.25]

equationNumber:[1056; 23]

Node

-coordinates : double[]

-dofsValues : double[]

-equationNumbers : int[]

-firstOrderTimeDerivatives:double[]

-secondOrderTimeDerivatives:double[]

-incrementsValues : double[]

-numbered : boolean

-valuesPredictors : double[]

+coordinates(double[] toto)

+coordinates()

+extractValuesFrom(MathVector V)

+…

Figure 113 Exemple de classe

1.2.2 Héritage

Un concept essentiel en programmation orienté objet est celui d‟héritage. La conception de la

nouvelle classe, qui hérite des propriétés et des aptitudes de l‟ancienne, peut ainsi s‟appuyer

sur des réalisations antérieures parfaitement au point et les spécialiser à volonté. L‟héritage

facilite la réutilisation de produits existants, d‟autant plus qu‟il peut être réitéré autant de fois

que nécessaire.

MVComponent

ImposedValues

BoundaryCondition

IncrementalBoundaryCondition

InitialCondition

SchwarzBoundaryCondition

Load

NodalLoad

SurfaceLoad

Figure 114 Exemple de hiérarchie de classes

Annexe 1

150

Dans la Figure 114, un exemple d‟une hiérarchie partielle de code élément fini est donné. Les

Classes ImposedValue et Load héritent de la superclasse MVComponent (Model View

Component) le comportement de gérer l'impression dans l'environnement graphique.

Le comportement particulier lié à la gestion des degrés de liberté est implémenté dans les

sous-classes de la classe ImposedValue. La gestion des charges (charges nodales ou charges

surfacique) est implémentée dans les sous-classes de la classe Load.

1.2.3 Polymorphisme

Le polymorphisme vient compléter la notion d‟héritage, il permet de manipuler des objets

sans en connaître complètement le type en se basant sur la relation d‟héritage. Le

polymorphisme se traduit par la compatibilité entre un type dérivé et le type de base. A cela

s‟ajoute la gestion du lien dynamique pour les méthodes, dans le sens où il permet d‟obtenir le

comportement adapté à chaque type d‟objet sans avoir besoin de tester sa nature. Un objet a

comme type non seulement sa classe mais aussi n‟importe quelle classe dérivée.

1.2.4 Interface

Java autorise la séparation entre la définition du comportement d'un objet et l‟implantation

correspondante. L'écriture d'une interface, puis d'une classe implémentant cette interface

réalise cette opération.

Une classe peut implémenter plusieurs interfaces alors qu‟une classe ne peut dériver que

d‟une seule classe abstraite. Une interface est une collection de définitions de méthodes (sans

implémentation) et de valeurs constantes. On utilise des interfaces pour définir un

comportement qui pourra ensuite implémenter n'importe quelle classe. Définir une interface

est similaire à définir une classe. Un exemple d‟interface est donné dans la Figure 115.

static public interface RadialReturnIntegrable extends

Integrator.Integrable {

public void computePlasticCorrection(Tensor

trialStress, GaussPoint gp, Element element);

}

Figure 115 Exemple d‟une déclaration d‟interface

Dans la Figure 116, la classe LinearElasticPerfectlyPlastic implémente l‟interface

RadialReturnIntegrable. La méthode de l‟interface est alors programmée dans la classe

LinearElasticPerfectlyPlastic.

Annexe 1

151

public class LinearElasticPerfectlyPlastic implements

RadialReturnAlgorithm.RadialReturnIntegrable {

…

public void computePlasticCorrection(Tensor trialStress,

GaussPoint gp, Element element) {

Tensor DfDsigma = criterion.firstOrderDerivative(trialStress);

FullMatrix Del = element.computeConstitutiveMatrix();

double k = material.getProperty("sigy") / Math.sqrt(3.0);

double f = criterion.value(trialStress, k);

double dGamma = f / (DfDsigma.times(Del.times(DfDsigma)));

Tensor dEpsp = DfDsigma.times(dGamma);

trialStress.minus(Del.times(dEpsp));

gp.plastic();

gp.getInternalParameters().get("dGamma").copyValue(dGamma);

}

…

}

Figure 116 Extrait d‟une classe implémentant une interface

1.2.5 Classes internes

Java autorise la définition d‟une classe à l‟intérieur d‟une autre classe. C‟est la notion de la

classe interne (inner class). Un des intérêts de cette notion est qu‟une classe interne peut avoir

accès aux méthodes et attributs de la classe englobante. L'accessibilité à la classe interne

respecte les règles de visibilité du langage. Il est même possible de définir une classe interne

privée pour limiter son accès à sa seule classe principale. Pour pouvoir utiliser une variable de

classe dans une classe interne, il faut la déclarer dans sa classe englobante. Il existe quatre

types de classes internes :

1.2.5.1 Les classes internes non statiques :

Les classes internes non statiques sont définies dans une classe dite "principale " en tant que

membre de cette classe. Leur avantage est de pouvoir accéder aux autres membres de la classe

principale même ceux déclarés comme privés. Dans le code, pour pouvoir faire référence à un

membre de la classe principale, il suffit simplement d'utiliser son nom de variable. Une classe

peut hériter d'une classe interne. Dans ce cas, il faut obligatoirement fournir aux constructeurs

Annexe 1

152

de la classe une référence sur la classe principale de la classe mère et appeler explicitement

dans le constructeur de cette classe principale.

1.2.5.2 Les classes internes statiques

Les classes internes statiques sont des classes internes qui ne possèdent pas de référence vers

leur classe principale. Elles ne peuvent donc pas accéder aux membres d'instance de leur

classe englobante. Elles peuvent toutefois avoir accès aux variables statiques de la classe

englobante. Leur utilisation est obligatoire si la classe est utilisée dans une méthode statique

qui par définition peut être appelée sans avoir d'instance de la classe et que l'on ne peut pas

avoir une instance de la classe englobante. Comme elle ne possède pas de référence sur sa

classe englobante, une classe interne statique est en fait traduite par le compilateur comme

une classe principale.

1.2.5.3 Les classes internes locales

Les classes internes locales sont définies à l'intérieure d'une méthode ou d'un bloc de code.

Ces classes ne sont utilisables que dans le bloc où elles sont définies. Les classes internes

locales ont toujours accès aux membres de la classe englobante. Leur particularité, en plus

d'avoir un accès aux membres de la classe principale, est d'avoir aussi un accès à certaines

variables locales du bloc où est définie la classe interne.

Dans la Figure 117, deux classes décrivant une formulation pour un problème découlement

(flux de Darcy) sont présentées. Sur le côté droit de la Figure 117, la classe

PressureDarcyPenalizedFormulation représente la formulation éléments finis. Ses attributs

correspondent aux différents paramètres d‟une formulation variationnelle : paramètre de

pénalisation, de la tolérance sur la position de surface libre, …. La méthode initialize() permet

de définir un champ discrétisé sur le domaine correspondant à l‟inconnue du problème. La

méthode getElement() permet de récupérer les éléments correspondant à la formulation.

La classe qui définit le problème est la classe PressureDarcy, elle-même est définie dans la

classe PressureDarcyPenalizedFormulation (sur le côté gauche de la Figure 117). La classe

PressureDarcy définit le calcul des contributions élémentaires nécessaires à l‟algorithme. Les

instances de cette classe peuvent avoir accès localement aux données de la formulation définie

dans la classe PressureDarcyPenalizedFormulation.

Annexe 1

153

Figure 117 Classe interne pour gérer une formulation éléments finis

1.2.5.4 Les classes internes anonymes

Les classes internes anonymes sont des classes internes qui ne possèdent pas de nom. Elles ne

peuvent donc être instanciées qu'à l'endroit où elles sont définies. Ce type de classe est très

pratique lorsqu'une classe doit être utilisée une seule fois. Les classes anonymes sont un

moyen pratique de déclarer un objet sans lui avoir trouvé un nom. La contrepartie est que

cette classe ne pourra être instanciée dans le code qu'à l'endroit où elle est définie : elle est

déclarée et instanciée en un seul et unique endroit.

1.2.6 Les paquetages (ou packages)

Les paquetages servent à structurer l'ensemble des classes et interfaces écrites en Java. Cette

notion est proche de celle de bibliothèque que l‟on rencontre dans d‟autres langages. Un

paquetage est caractérisé par un nom qui est soit un simple identificateur, soit une suite

d‟identificateurs séparés par des points.

La hiérarchie d'un paquetage se retrouve dans l'arborescence des répertoires physiques.

Chaque paquetage est dans un répertoire nommé du nom du paquetage. D'une façon générale,

l'instruction package associe toutes les classes qui sont définies dans un fichier source à un

même package. Il existe plusieurs types de paquetages : le paquetage par défaut (identifié par

le point qui représente le répertoire courant et permet de localiser les classes qui ne sont pas

associées à un paquetage particulier), les paquetages standards qui sont empaquetés dans le

Annexe 1

154

fichier classes.zip et les paquetages personnels. Un paquetage par défaut est systématiquement

attribué par le compilateur aux classes qui sont définies sans déclarer explicitement une

appartenance à un paquetage. Ce paquetage par défaut correspond au répertoire courant qui

est le répertoire de travail.

Dans la Figure 118 La classe Matrix, FullMatrix et Vector appartiennent au package

fem.linearalgebra, et la classe MisesPerfectlyPlastic appartient au package fem.material.

Dans cet exemple, la classe FullMatrix peut accéder directement à certains attributs de la

classe Vector afin d'améliorer l'efficacité de l‟implémentation en évitant les appels

supplémentaires des méthodes, par exemple dans le produit entre matrice et vecteur. Ils

appartiennent à un même package, et un accès partiel a été donnée dans ce contexte. En

dehors du package linearalgebra, ce partage partielle n'est pas autorisée, par exemple la

classe MisesPerfectlyPlastic n'a pas accès aux éléments de la matrice. L'organisation des

packages suivent l'organisation des répertoires. Ainsi, l‟implémentation du code suit

l'organisation abstraite.

fem

linearalgebra

Matrix

FullMatrix

Vector

material

MisesPerfectlyPlastic

…

Légende:

-Les packages sont représentés en italique

-Les classes sont représentées en écriture gras.

Figure 118 Une hiérarchie de packages et les classes associées

1.3 Des paquetages particuliers pour la simulation numérique en

mécanique

Au delà des aspects langage de programmation, Java est un environnement conçu pour

appréhender toutes les problématiques de l‟informatique moderne : communication réseau,

distribution, gestion réseau, interfaçage home machine, procédures graphique 3D,…

Annexe 1

155

1.3.1 RMI (Remote Method Invocation)

RMI est une API Java permettant de manipuler des objets distants (c'est-à-dire un objet

instancié sur une autre machine virtuelle, éventuellement sur une autre machine du réseau) de

manière transparente pour l'utilisateur, c'est-à-dire de la même façon que si l'objet était sur la

machine virtuelle (JVM) de la machine locale. Ainsi, un serveur permet à un client d'invoquer

des méthodes à distance sur un objet qu'il instancie. Deux machines virtuelles sont donc

nécessaires (une sur le serveur et une sur le client) et l'ensemble des communications se fait

en Java.

Figure 119 Un exemple (dans FEMJava) de diagramme de classe

La Figure 119 montre que l'ordre des calculs est globalement le même que dans la version

multithread. La seule différence consiste dans la couche supplémentaire en charge des cadres

simultanés. L'idée principale est que les sous-domaines, qui sont censés être exécutés

simultanément, sont incorporés dans un objet simultané. Le cadre du travail RMI permet à un

objet d‟invoquer des méthodes d'objets distants à partir d'autres objets situés dans d'autres

machines virtuelles Java.

1.3.2 Gestion des processus légers

Certaines applications, y compris scientifiques, nécessitent un haut degré de parallélisme. Les

threads permettent d‟exploiter efficacement les machines multiprocesseurs et la répartition

Annexe 1

156

des tâches. Le développement à l‟aide de threads peut apporter plusieurs exécutions

simultanées des tâches. On ne peut pas avoir plusieurs threads qui essayent d‟accéder au

même variable sans coordination. Java rend la création, le contrôle et la coordination des

threads très simple, en intégrant quelques uns de ces concepts dans le langage.

CODE A------------------------------------------------------ Thread thread = new Thread ( new Runnable () { public void run()

{ // Cette partie du programme doit être executer // simultanément dans l’application }} ) ;

thread.start(); // Exécuter le THREAD

CODE B------------------------------------------------------ for (int i = 0; i < size; i++) {

for (int j = 0; j < coefficients[i].length; j++) {

Vanswer[i]+=(coefficients[i][j]*vdouble[profile[i][j]-1]);

}}

CODE C------------------------------------------------------

for ( int ii = 0 ; ii < numberOfThreads ; ii++ ){

final int Ni = ii ;

Thread threadii = new Thread( new Runnable(){

int number = Ni ;

public void run(){

for(int i = bounds[number][0];i<bounds[number][1]; i++){

for ( int j = 0 ; j < coefficients[i].length ; j++ ){

Vanswer[i] += ( coefficients[i][j] *

vdouble[profile[i][j] - 1] ) ;}}}}) ;

productThreads[ii] = threadii ; }

for ( int i = 0 ; i < numberOfThreads ; i++ )

productThreads[i].start() ; // Commencer le THREAD i

Figure 120 Différentes implémentations du produit Matrice/Vecteur

Dans la Figure 120, le Code A est une implémentation typique d'un objet Thread. Cette

implémentation a l'avantage d'être exécuté en même temps que le code et d'introduire la

méthode run () dans une classe anonyme qui implémente l‟interface Runnable. Dans le code

Annexe 1

157

B et C, une comparaison est faite entre un simple produit matrice/vecteur dans le cas d'une

matrice creuse.

1.3.3 Java 3D et interfaces graphiques (GUI)

Dans les programmes à interfaces graphiques, la communication avec l‟utilisateur se fait par

l‟intermédiaire de composants tels que les menus déroulants, les menus surgissant, les barres

d‟outils ou les boîtes de dialogues, ces dernières pouvant renfermer des composants aussi

variés que les boutons poussoirs, les cases à cocher, les boutons radio, les boîtes de saisie,…

Des paquetages en Java contiennent ce type de composant pour construire des interfaces

graphiques ergonomiques.

Figure 121 Exemple d‟application graphique intégrant du graphisme 3D

Java 3D est un API Java destiné à l'affichage 3D. Java 3D est destiné à l'écriture

d'applications et d'Applets graphiques 3D. Ses principales fonctionnalités sont : la

modélisation et l‟affichage 3D de scènes, l‟évolutivité par dérivation de classes, la possibilité

d'implanter des comportements spécifiques programmés, la détection des collisions entre

objets,… (voir Figure 121).

Annexe 1

158

Annexe 2

159

Annexe 2 Dérivation des équations de

Navier-Stokes en milieu incompressible

2.1 Problème d’écoulement d’un fluide

2.1.1 Les équations d’équilibre

On considère la formulation de Navier-Stokes pour un fluide newtonien incompressible. On

cherche la vitesse et la pression telles que les équations d‟équilibre, les conditions de

bords et les conditions d‟incompressibilité soient respectées.

Les conditions de bords décrites dans la Figure 122.

Figure 122 Définition du domaine et des conditions de bord

2

1

f

F

uu

Annexe 2

160

Les équations du problème sont :

2.1.2 Formulation variationnelle

On se place en régime permanent, la formulation éléments finis utilisée pour développer le

modèle numérique est une formulation mixte (vitesse-pression), formulation de Galerkin

stabilisée par ajout de termes de type moindres-carrés (voir [Tezduyar 92]. Elle est donnée

comme suit :

Considérons les ensembles suivants :

La formulation variationnelle approximée du problème est :

où le paramètre de stabilisation est de la forme :

2.2 Génération des contributions élémentaires

On peut exprimer cette formulation en fonction du formalisme développé dans le Chapitre 3.

On obtient les huits termes intégraux suivants :

Annexe 2

161

Le terme (terme de convection) conduit à écrire pour la contribution élémentaire un

tenseur d‟ordre 2 de la forme :

En identifiant les indices:

, , , ,

,

,

,

,

;

;

Donc:

Annexe 2

162

Le tenseur est de dimension est de la forme

On effectue le même développement pour les termes jusqu‟au pour obtenir

respectivement les tenseurs jusqu‟au

.

,…,

sont les contributions élémentaires de la formulation mixte proposée.

2.3 Description symbolique du problème

Dans cette partie, l‟outil de discrétisation automatique est appliqué sur les champs de la

formulation afin de construire les contributions élémentaires.

Dans les Figure 123, on donne le script permettant de construire les éléments de base de la

formulation.

Annexe 2

163

Field u = new Field("u", new int[] { 3 }, 3, true, new

Operator[] {new Field("N1", null), new Field("N2", null),



new Field("N7", null), new Field("N8", null) });

DifferentialOperator gradu=new DifferentialOperator("grad",u);

DifferentialOperator Eu = new DifferentialOperator("E", u);

DifferentialOperator divu=new DifferentialOperator("div", u);

Field p = new Field("p", null, 3, true, new Operator[] {




new Field("N7", null), new Field("N8", null) });

DifferentialOperator gradp=new DifferentialOperator("grad",p);

Field w = new Field("w", new int[] { 3 }, 3, true,null);

DifferentialOperator divw = new DifferentialOperator("div",w);

DifferentialOperator gradw=new DifferentialOperator("grad",w);

DifferentialOperator Ew = new DifferentialOperator("E", w);

Field q = new Field("q", null, 3, true,null);

DifferentialOperator gradq=new DifferentialOperator("grad",q);

Field f = new Field("f", new int[] { 3 }, 3);

Fraction rho = new Fraction(100.0,1.0) ;

Fraction mu = new Fraction(1.0,1.0) ;

Tensor Tq = q.giveDiscretizationTensor();

Tensor Tp = p.giveDiscretizationTensor();

Tensor Tw = w.giveDiscretizationTensor();

Tensor Tgradu = gradu.giveDiscretizationTensor();

Tensor Tgradw = gradw.giveDiscretizationTensor();

Tensor Tdivu = divu.giveDiscretizationTensor();

Tensor Tdivw = divw.giveDiscretizationTensor();

Figure 123 La construction du problème

La construction des contributions élémentaires est montrée dans les Figure 124 et Figure 125

//-----------------------------A1-----------------------------

Tensor TuTgradu = Tu.simpleContraction(Tgradu);

Tensor TugraduTw = TuTgradu.simpleContraction(Tw);

Tensor rhoTugraduTw = rho.scalarProduct(TugraduTw);

//-----------------------------A2-----------------------------

Tensor TEu = Eu.giveDiscretizationTensor();

Tensor TEw = Ew.giveDiscretizationTensor();

Tensor TEuTEw = TEu.doubleContraction(TEw);

Fraction deux = new Fraction(2,1);

Fraction deuxmu = deux.multiplication(mu);

Tensor deuxmuTEwTEu=deuxmu.scalarProduct(TEuTEw);

Figure 124 Contributions élémentaires pour le problème de Navier-Stoke (partie 1)

Annexe 2

164

//-----------------------------A3----------------------------

Tensor TdivuTq = Tdivu.tensorProduct(Tq);

//-----------------------------A8----------------------------

Tensor tauxTgradqTgradp = tauxTgradq.simpleContraction(Tgradp);

Tensor unSurRhotauxTgradqTgradp =

tauxTgradqTgradp.scalarProduct(unSurRho);

//----------------------Tenseurs------------------------

Tensor[] TSymbolicLeft = new Tensor[] { rhoTugraduTw,

deuxmuTEwTEu, TdivuTq, moinspTdivw};

Tensor[] TSymbolicright = new Tensor[] { moinsfv};

//----------------------Formulation------------------------

Formulation g = new Formulation("NavierStockes", new

Equation[]{equation}, new Material (lambda, mu),


Figure 125 Contributions élémentaires pour le problème de Navier-Stoke (partie 2)

2.3 Application numérique

On considère le problème de la cavité entraînée en 3D donné dans la Figure 126. Le domaine

est cubique d‟arête . On impose une vitesse nulle sur toutes les faces sauf sur la dernière

pour laquelle une vitesse parallèle à une des arêtes d‟intensité . Les données du

matériau sont : la densité volumique et le module de cisaillement .

Le nombre de Reynolds pour le problème est donc .

Figure 126 Cavité entraînée en 3D

Dans la Figure 127, on montre la solution numérique du problème en vitesse et en pression.

v=1

v=0

v=0 v=0

Annexe 2

165

Figure 127 Solution en vitesse et en pression

Annexe 2

166

Références

167

Références

[Adeli 00] Adeli H, Kim H, “Web-based interactive courseware for structural steel design

using Java” Computer Aided Civil and Infrastructure Engineering, 15 (2), 158-166, 2000.

[Baduel 04] Baduel L., F. Baude, D. Caromel, C. Delbé, N. Gama, S. El Kasmi and S.

Lanteri, “A parallel object-oriented application for 3-D electromagnetism”, ECCOMAS,

Jyväskylä, Finland 2004.

[Barbier 92] Barbier C., “Automatic generation of bending element matrices for finite element

method using REDUCE”, Engineering Computations, 9, 477-494, 1992.

[Bathe 96] Bathe K-J, “Finite element procedure”, Simon and Schuster/A Viacom Company,

Upper Saddle River, New Jersey, 1996.

[Baugh 92] Baugh J. W., Rehak D. R., “Data Abstraction in Engineering Software

Development”, Journal of Computing in Civil Engineering, 6 (3), July, 1992.

[Besson 97] Besson J. and R. Foerch, “Large scale object-oriented finite element code

design”, Comput. Methods Appl. Mech. Engrg., 142, 165-187, 1997.

[Bull 01] Bull J.M., L. A. Schmith, L. Pottage and R. Freeman, “Benchmarking Java against

C and Fortran for Scientific Applications”, Joint ACM JavaGrande – ISCOPE 2001

Conference, Stanford Universtity, June 2-4, 2001.

[Cameron 97] Cameron F., “Automatic generation of efficient routines for evaluating

multivariate polynomials arising in finite element computations”, Adv. in Engr. Soft., 28,

239-245, 1997.

[Cecchi 77] Cecchi M.M. and C. Lami, “Automatic generation of stiffness matrices for finite

element analysis”, Internat. J. Numer. Methods Engrg., 11, 396-400, 1977.

[Choi 92] Choi DK., Nomura, S. “Application of symbolic computation to two-dimensional

elasticity”, Computers & Structures, 43, 645-649, 1992.

Références

168

[Devloo 94] Devloo P.R.B., “Efficient issues in an object oriented programming

environment”, Proceedings of CST 94 , Athens, Greece, vol. Artificial intelligence and object

oriented approaches for structural engineering, Civil Comp Press, 147-151, 1994.

[Dubois-Pelerin 93] Dubois-Pelerin Y. and Th. Zimmermann, “Object-oriented finite element

programming: III – An efficient implementation in C++”, Computer Methods Appl. Mech.

Eng., 108(2), 165-183, 1993.

[Eyheramendy 96] Eyheramendy D., Th. Zimmermann, “The Object-oriented finite elements:

II. A symbolic environment for automatic programming”, Comput. Methods Appl. Mech.

Engrg., 32, 259-276, 1996.

[Eyheramendy 98] Eyheramendy D., Th. Zimmermann, “Object-oriented finite elements: III.

Theory and application of automatic programming”, Comput. Methods Appl. Mech. Engrg.,

154, 41- 68, 1998.

[Eyheramendy 03] Eyheramendy D., “Object-Oriented parallel CFD with Java”, Moscow

Russia, Eds. Chetverushkin, Ecer, Satofuka, Périaux, Fox, Ed. Elsevier, 409-416, 2003.

[Eyheramendy 04] Eyheramendy D. and D. Guibert, “A Java Approach for Finite Elements

Computational Mechanics”, ECCOMAS 2004, Jyvaskyla, Finland, July 2004.

[Eyheramendy 06] Eyheramendy D., “High abstraction level frameworks for the next decade

in computational mechanics”, Innovation in Engineering Computational Technology, Eds.

B.H.V. Topping, G. Montero and R. Montenegro, ©Saxe-Cobourg Publications, Chap. 3, 41-

61, 2006.

[Eyheramendy 07] Eyheramendy D. and F. Oudin, “Advanced object-oriented techniques for

coupled multiphysics”, In Civil Engineering Computation: Tools and Techniques, Ed. B.H.V.

Topping, ©Saxe-Cobourg Publications, ISBN 978-1-874672-32-6, Chap. 3, 37-60, 2007.

[Eyheramendy 08] Eyheramendy D. and F. Oudin, “Object-oriented finite elements: From

Smalltalk to Java”, In Trends in Engineering Computational Technology Eds. M.

Papadrakakis and B.H.V. Topping, ©Saxe-Cobourg Publications, 2, 17-40, 2008.

[Eyheramendy 11] Eyheramendy D., R. Saad et F. Oudin “Advanced object-oriented

paradigms for parallel computational mechanics”, Proceedings of the 2nd

International

Conference on Parallel, Distributed, Grid and Cloud Computing for Engineering, B.H.V.

Références

169

Topping, J.M. Adam, F.J. Pallarés, R. Bru and M.L. Romero, (Editors), Civil-Comp Press,

Stirlingshire, Scotland, 2011.

[Fang 02] Fang Y., et al., “Influence of surface residual stress state on crack path evolution in

polycrystalline alumina”, J Am Ceram, 85(7), 1783–1787, 2002.

[Fenves 79] Fenves S.J. and Korncoff A.R, “Symbolic generation of finite element stiffness

matrices”, Computers & Structures, 10, 119-124, 1979.

[Fenves 90] Fenves G.L., “Object-Oriented programming for engineering software

development”, Engineering with Computer, 6, 1-15, 1990.

[Flanagan 04] Flanagan D., Java Examples in a Nutshell, Ed. O'Reilly, 2004.

[Flanagan 06] Flanagan D, Java in a Nutschell, Ed. O'Reilly, 2006.

[Foerch 95] Foerch R., J. Besson, P. Pilvin and G. Cailletaud, “Formulation des relations de

comportement dans les calculs par éléments finis : approche C++”, Actes du 2nd Colloque

national en calcul des structures, Giens, Hermès, 547-552, 1995.

[Forde 90] Forde B.W.R, R.O. Foschi and S.F Stiemer, “Object-Oriented Finite Element

Analysis”, Computers & Structures, 34, 355-374, 1990.

[Fritzson 92] Fritzson, P., Fritzson D., “The need for highlevel programming support in

scientific computing applied to mechanical analysis”, Computers & Structures, 45, 387-395,

1992.

[Garrigues 07] Garrigues J., “Fondements de la mécanique des milieux continus”,

éditeur : Hermes Sciences, 2007.

[Golay 02] Golay F., Breitkopf P., Touzot G., “SIC :Système Interactif de Conception”,

éditeur :Université de Technologie de Compiègne, 2002.

[Gunderson 71] Gunderson R.H. and A. Cetiner, “Element stiffness matrix generator”, J.

Struct. Div., Poceedings of ASCE, 363-375, 1971.

[Häuser 00] Häuser J., T. Ludewig, R.D. Williams, R. Winkelmann, T. Gollnick, S. Brunett

and J. Muylaert, “A test suite for high-performance parallel Java”, Advances in Engineering

Software, 31, 687-696, 2000.

Références

170

[Heng 06] Heng B.C.P. and R.I. Mackie, “Design Patterns in Object-Oriented Finite Element

Programming”, in Proceedings of the Fifth International Conference on Engineering

Computational Technology, B.H.V. Topping, G. Montero and R. Montenegro, (Editors),

Civil-Comp Press, United Kingdom, 2006.

[Hoa 80] Hoa S.V. and S. Sankar, “A computer program for automatic generation of stiffness

and mass matrices in finite-element analyis”, Computers & Structures, 11, 147-161, 1980.

[Hughes 78] Hughes T.J.R, K.S Pister, “Consistent linearization in mechanics of solids and

structures”, Comput. Struct., 8, 391-397, 1978.

[Ibrahimbegovic 09] Ibrahimbegovic A., “Nonlinear solid mechanics: theoretical

formulations and finite element solution methods”, Springer, Berlin, 1-571, 2009.

[Ioakimidis 93] Ioakimidis NI. “Elementary applications of MATHEMETICA to the solution

of elasticity problems by the finite element method”, Comput. Methods Appl. Mech. Engrg.,

102, 29-40, 1993.

[Korelc 97] Korelc J., “Automatic generation of finite-element code by simultaneous

optimization of expressions”, Theoretical Computer Science, 187, 23 l-248, 1997.

[Korelc 02] Korelc J., “Multi-language and Multi-environment Generation of Nonlinear

Finite Element Codes”, Engineering with Computers, 18, 312-327, 2002.

[Korelc 09] Korelc J., “Automation of primal and sensitivity analysis of transient coupled

problems”, Compt. Mech., 44, 631-649, 2009.

[Leff 91] Leff L. and Y.Y. Yun, “The symbolic finite element analysis system”, Computers &

Structures, 41, 227-231, 1991.

[Liu 01] Liu YQ, “Fast solution for large scale linear algebraic equations in finite element

analysis”, Acta Mechanica Solida Sinica, 14, 89-94, 2001.

[Logg 07] Logg A., “Automating the Finite Element Method”, Arch Comput Methods Eng.,

14, 93-138, 2007.

[Luft 71] Luft W., J.M. Roesset and J.J. Connor, “Automatic generation of finite element

matrices”, Struct. Div. Proceedings of ASCE, 349-361, 1971.

Références

171

[Ma 02] Ma YQ., “Object-oriented finite element analysis and programming in VC++”, Appl

Math Mech, 23(12), 1437–1443, 2002.

[Mackerle 04] Mackerle J., “Object-oriented programming in FEM and BEM: a bibliography

(1990-2003)”, Advances in Engineering Software, 35, 325-336, 2004.

[Mackie 00] Mackie RI., “An object-oriented approach to calculation control in finite element

programs”, Comput Struct, 77(5), 461–474, 2000.

[Mackie 01] Mackie RI., “Implementation of sub-structuring within an object-oriented

framework”. Adv Eng Software, 32(10), 749–758, 2001.

[Mackie 01] Mackie RI., “Using objects to handle calculation control in finite element

analysis”, Dev Eng Comput Tech, Edinburgh: Civil-Comp, 123–130, 2000.

[Mackie 07] Mackie R.I., “Object-oriented design of pre-conditionned iterative equation

solvers using .NET”, Proceedings of 5th Int. Conf. on Engineering Computational

Technology, Las Palmas de Gran Canaria, Spain, 12-15 Sept. 2007.

[Menétrey 93] Menétrey Ph. and Th. Zimmermann, “Object-Oriented Non-Linear Finite

Element Analysis : Application to J2 plasticity”, Computers & Structures, 49 , 5, 767-777,

1993.

[Miller 98] Miller G.R, “A LISP-Based Object-Oriented approach to structural analysis”,

Engineering with Computer, 4, 197-203, 1988.

[Nikishkov 99] Nikishkov G.P. and H. Kanda, “The development of a Java engineering

application for higher-order asymptotic analysis of crack-tip fields”, Advances in

Engineering Software, 30, 469-477, 1999.

[Nikishkov 03] Nikishkov G.P., “Generating contours on FEM/BEM higher-order surfaces

using Java 3D textures”, Advances in Engineering Software, 34, 469-476, 2003.

[Nikishkov 06] Nikishkov G.P., “Object-oriented design of a finite element code in Java”,

Computer Modeling in Engng and Sciences, 11, 81-90, 2006.

[Nikishkov 10] Nikishkov G.P., “Programming Finite Elements in Java”, Springer Eds.

ISBN, 978-1-84882-971-8, 2010.

Références

172

[Noor 81] Noor A.K. and C.M. Andersen, “Computerized symbolic manipulation in nonlinear

finite element analysis”, Computers & Structures, 13, 379-403, 1981.

[Padial-Collins 04] Padial-Collins N.T., W.B. VanderHeyden, D.Z. Zhang, E.D. Dendy and

D. Livescu, “Parallel operation of CartaBlanca on shared and distributed memory

computers”, Concurrency and Computation: Practice and Experience 16, 61-77, 2004.

[Rehak 89] Rehak D.R. and J.W. Baugh Jr., “Alternative Programming Techniques for Finite

Element Programming Development”, Proceedings IABSE Colloquium on Expert Systems in

Civil Engineerings, Bergamo, Italy. IABSE, 1989.

[Riley 03] Riley C.J., S. Chatterjee and R. Biswas, “High-performance Java codes for

computational fluid dynamics”, Concurrency and Computation: Practice and Experience 15,

395-415, 2003.

[Rio 08] Rio G., Laurent H., Bles G., “Asynchronous interface between a finite element

commercial software ABAQUS and an academic research code HEREZH”, Adv. in Engr.

Soft, 39, 1010-1022, 2008.

[Scholz 92] Scholz S. P., “Elements of an Object-Oriented FEM++ Program in C++”,

Computers & Structures, Vol. 43, No. 3, 517-529, 1992.

[Silvester 94] Silvester, PP., Chamlian SV., “Symbolic generation of finite elements for skin-

effect integral equations”. IEEE Trans. Magnetics, 30(5), 3594-3597, 1994.

[Tezduyar 92] Tezduyar TE., “finite element formulations for incompressible flow

computations”, Advances in applied Mechanics, 28, 1-44, 1992.

[VanderHeyden 03] VanderHeyden W.B., E.D. Dendy and N.T. Padial-Collins,

“CartaBlanca-a pure-Java, component-based systems simulation tool for coupled nonlinear

physics on unstructured grids-an update”, Concurrency and Computation: Practice and

Experience, 15, 431-458, 2003.

[Verpeaux 88] Verpeaux P., T. Charras et A. Millard, “CASTEM 2000 : une approche

moderne du calcul de structures”, Calcul des structures et intelligence artificielle, Vol. 2, J.M.

Fouet, P. Ladevèze et R. Oyahon Eds., Pluralis, 1988.

Références

173

[Yagawa 90] Yagawa, G; Ye, GW; Yoshimura, S, “A numerical integration scheme for finite

element method based on symbolic manipulation”, Int. J. Numer. Methods Engrg., 29, 1539-

1549, 1990.

[Yang 94] Yang, CY. “An algebraic-expressed finite element model for symbolic

computation”, Computers & Structures, 52(5), 1069-1077, 1994.

[Zeglinski 94] Zeglinski G. W., Han R. P. S.; Aitchison P., “Object-Oriented Matrix Classes

For Use In A Finite Element Code Using C++”, International Journal For Numerical Methods

In Engineering, 30, 3921-3937, 1994.

[Zimmerman 92] Zimmermann Th., Y. Dubois-Pèlerin and P. Bomme, “Object-oriented finite

element programming: I. Governing principles”, Comput. Methods Appl. Mech. Engrg., 98,

291-303, 1992.

[Zimmerman 96] Zimmerman Th. And Bomme P, “Towards intelligent objects in finite

element programming”, Proc Third Int Conf on Computational Structures Technology,

Budapest, 107-114, 1996.

Documents

sur une approche à objets généralisée pour la mécanique non linéaire