RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

SUR L'IMPÉLÉMENTATION
DE PROBLÈMES GÉNÉRAUX
DE COQUES

Michel BERNADOU
Jean-Marie BOISSERIE

Rapport de Recherche N° 317

Juillet 1978
SUR L'IMPLEMENTATION DE PROBLEMES GENERAUX DE COQUES

M. Bernadou - J.M. Boisserie

RÉSUMÉ

Dans ce rapport, nous détaillons l'implémentation du modèle linéaire général de coques de W.T. Koiter. La surface moyenne de la coque est définie à l'aide d'une carte utilisant des coordonnées curvilignes quelconques. Pour approcher la solution, nous utilisons une méthode conforme d'éléments finis ainsi que des schémas d'intégration numérique. Ce travail est complété par la résolution numérique de problèmes particulièrement simples dont la solution exacte est connue. Les listes FORTRAN correspondantes sont données en annexe.

ABSTRACT

In this report, the implementation of Koiter's general linear shell model is considered. The middle surface of the shell is defined by a mapping using general curvilinear coordinates. To approximate the solution, a conforming finite element method is used in conjunction with numerical quadrature schemes. This work is completed by numerical examples of particularly simple problems for which the exact solution is known. The associated FORTRAN listing is given in the appendix.
# INTRODUCTION

Après avoir souligné qu'un des aspects extrêmement intéressant de la méthode des éléments finis est son caractère modulaire, nous avons développé dans un travail précédent (cf. BERNADOU-BOISSERIE [1]) les modules suivants :

(i) module d'interpolation à l'aide de l'élément fini d'ARGYRIS. En particulier, nous montrons comment l'utilisation conjointe des coordonnées barycentriques et des paramètres d'excentricité permet de définir l'élément fini pour un triangle quelconque, indépendamment de toute utilisation d'un élément fini de référence ;

(ii) module "fonctionnelle d'énergie" pour un problème de potentiel ;

(iii) module "fonctionnelle d'énergie" pour un problème de plaque. Pour ce problème, le module d'interpolation ci-dessus, réalisé à l'aide de l'élément fini d'ARGYRIS, trouve sa pleine utilisation : il permet en effet de construire un sous-espace de dimension finie de l'espace de Sobolev $H^2(\Omega)$, i.e., d'obtenir une méthode conforme d'éléments finis.

Comme suite à ce travail, nous proposons dans ce rapport un module "fonctionnelle d'énergie de coques" correspondant au modèle général de coques de W.T. KOITER [1,2]. Ce modèle est formulé sur la surface moyenne définie à l'aide d'une carte utilisant des coordonnées curvilignes quelconques. Pour réaliser ce module nous utilisons le module d'interpolation mentionné ci-dessus ainsi qu'un schéma d'intégration numérique satisfaisant les critères établis dans BERNADOU [1,2].

Dans la Section 3.1 nous indiquons comment calculer la matrice de rigidité élémentaire lorsque les trois composantes du déplacement sont approchées dans un même espace d'éléments finis construit à l'aide de l'élément d'ARGYRIS. Dans la Remarque 3.1-1, nous précisons les modifications qu'il convient d'apporter à ce calcul lorsque les composantes tangentielles $u_1$ et $u_2$, d'une part, la composante normale $u_3$, d'autre part, sont approchées dans des espaces d'éléments finis différents. Puis, dans la Section 3.2 nous résolvons numériquement des problèmes particulièrement simples dont on connaît la solution exacte.
## Notations

La numérotation suit directement celle de BERNADOU-BOISSERIE [1] où se trouvent les deux premiers paragraphes. Par ailleurs les renvois à BERNADOU [1] sont indiqués par (II,...).

## Plan

3. Résolution numérique de problèmes de coques

3.1. Calcul de la matrice de rigidité et du second membre

3.2. Résolution numérique d'un problème particulier.

Annexe A

Annexe B

Annexe C

C.1. Module géométrie différentielle

C.2. Module rigidité élémentaire

C.3. Module second membre élémentaire
3 - RESOLUTION NUMERIQUE DE PROBLEMES DE COQUES

3-1. Calcul de la matrice de rigidité et du second membre :

Suivant (II, (3.1-1) (3.1-11) (3.1-12)) le problème discret associé
à la prise en compte de l'intégration numérique s'énonce :

$$
\text{(3.1-1)} \quad \left\{
\begin{array}{l}
\text{Trouver } \vec{u}_h \in \vec{V}_h \text{ tel que pour tout } \vec{v}_h \in \vec{V}_h, \\
\text{a}_h(\vec{u}_h, \vec{v}_h) = f_h(\vec{v}_h),
\end{array}
\right.
$$

où

$$
(3.1-2) \quad a_h(\vec{u}_h, \vec{v}_h) = \sum_{K \in \mathcal{K}} \sum_{\ell=1}^{L} \omega_{\ell,K} \, {}^t U_h(b_{\ell,K}) [A_{IJ}(b_{\ell,K})] V_h(b_{\ell,K}),
$$

$$
(3.1-3) \quad f_h(\vec{v}_h) = \sum_{K \in \mathcal{K}} \sum_{\ell=1}^{L} \omega_{\ell,K} \, {}^t F(b_{\ell,K}) V_h(b_{\ell,K}).
$$

Par analogie avec la relation (1.4-9) nous désignons par

$$
\left\{
\begin{array}{l}
\mathrm{DLGL}(u_1) \quad \text{la matrice des degrés de liberté globaux relatifs à } u_1, \\
\mathrm{DLGL}(u_2) \quad \text{la matrice des degrés de liberté globaux relatifs à } u_2, \\
\mathrm{DLGL}(u_3) \quad \text{la matrice des degrés de liberté globaux relatifs à } u_3.
\end{array}
\right.
$$

Alors, la matrice $^t U_h$, i.e.,

$$
t[u_1h \ u_1h,x \ u_1h,y \ u_2h \ u_2h,x \ u_2h,y \ u_3h \ u_3h,x \ u_3h,y \ u_3h,xx \ u_3h,xy \ u_3h,yy]
$$

s'écrit :
$$
(3.1-4) \quad t_{U_h} = [DG(\vec{u}_h)] \quad 1 \times 63 \quad
\left[
  \begin{array}{ccc}
    DA & O & O \\
    O & DA & O \\
    O & O & DA
  \end{array}
\right]
\left[
  \begin{array}{c|c}
    \lambda \quad \partial_x \lambda \quad \partial_y \lambda & O \\
    O & \lambda \quad \partial_x \lambda \quad \partial_y \lambda & O \\
    O & O & \lambda \quad \partial_x \lambda \quad \partial_y \lambda \quad \partial_{xx} \lambda \quad \partial_{xy} \lambda \quad \partial_{yy} \lambda
  \end{array}
\right] \quad (63 \times 63) \quad (63 \times 12)
$$

où

$$
(3.1-5) \quad [DG(\vec{u}_h)]_{1 \times 63} = [DLGL(u_{1h});DLGL(u_{2h});DLGL(u_{3h})].
$$

Suivant (II, Corollaire 3.5-1 et Figure 3.6-1) il convient d'utiliser
un schéma d'intégration exact pour les polynômes de degré 8, par exemple le
schéma d'intégration indiqué sur la Figure 2.1-1. Dans ce schéma, défini sur
un triangle quelconque, les noeuds d'intégrations sont repérés par leurs
coordonnées barycentriques. Soit $\omega_\ell$, $\ell = 1, \dots, L$, les poids correspondants.
Posons

$$
(3.1-6) \quad [DT] = \begin{bmatrix} DA & O & O \\ O & DA & O \\ O & O & DA \end{bmatrix} \qquad 63 \times 63
$$

et

$$
(3.1-7) \quad [\text{LAMBD}] = \left[ \begin{array}{c|c|c} \lambda & \partial_x \lambda & \partial_y \lambda \\ \hline O & \lambda \partial_x \lambda \partial_y \lambda & O \\ O & O & \lambda \partial_x \lambda \partial_y \lambda \partial_{xx} \lambda \partial_{xy} \lambda \partial_{yy} \lambda \end{array} \right] . \qquad (63 \times 12)
$$

Alors, les relations (3.1-2) et (3.1-3) s'écrivent respectivement

$$
(3.1-8) \quad a_h(\vec{u}_h, \vec{v}_h) = \sum_{K \in \mathcal{K}_h} \left\{ \text{aire}(K) [DG(\vec{u}_h)] [DT] \right\} \left\{ \left( \sum_{\ell=1}^{16} \omega_\ell [\text{LAMBD}] [A_{ij}]^t [\text{LAMBD}] \right) (b_\ell) \right\}^t [DT]^t [DG(\vec{v}_h)] ,
$$

$$
(3.1-9) \quad f_h(\vec{v}_h) = \sum_{K \in \mathcal{K}_h} \left\{ \text{aire}(K) \left\{ \sum_{\ell=1}^{16} \omega_\ell \, {}^t F(b_\ell) \, {}^t [\text{LAMBD}] (b_\ell) \right\} \, {}^t [DT] \, {}^t [DG(\vec{v}_h)] \right\} .
$$

Remarque 3.1-1. Si l'on veut utiliser des éléments finis distincts pour ap-
procher u₁, u₂ d'une part, u₃ d'autre part, les modifications portent sur
les matrices [DG], [DT] et [LAMBD].
Par exemple, si pour approcher $u_1$, $u_2$ on utilise l'élément fini d'Hermite de type 3, on a alors par analogie avec (1.6-1) :

$$u_{ch} = [DLGL1(u_{ch})]_{1 \times 10} [D1]_{10 \times 10} [A1]_{10 \times 10} [\lambda1]_{10 \times 1}, \alpha = 1,2.$$

Les matrices $[DLGL1(u_{ch})]$, $[D1]$, $[A1]$, $[\lambda1]$, sont définies par analogie avec les matrices $[DLGL(v_h)]$, $[D]$, $[A]$, $[\lambda]$ considérées en (1.4-9)(1.4-15) et (1.5-1).

La composante $u_3$ étant toujours approchée à l'aide de l'élément d'ARGYRIS, on aura cette fois-ci :

$$[DG(\vec{u}_h)]_{1 \times 41} = [DLGL1(u_{1h}); DLGL1(u_{2h}); DLGL(u_{3h})]$$

$$[DT]_{41 \times 41} = \begin{bmatrix} [D1] [A1] & 0 & 0 \\ 0 & [D1] [A1] & 0 \\ 0 & 0 & [D] [A] \end{bmatrix}_{41 \times 41}$$

$$[LAMBD] = \begin{bmatrix} \lambda_1 \partial_x \lambda_1 \partial_y \lambda_1 & 0 & 0 \\ 0 & \lambda_1 \partial_x \lambda_1 \partial_y \lambda_1 & 0 \\ 0 & 0 & \lambda_1 \partial_x \lambda_1 \partial_y \lambda_1 \partial_{xx} \partial_{xy} \partial_{yy} \lambda_1 \end{bmatrix}_{41 \times 12}$$

### 3-2. Résolution numérique d'un problème particulier

Dans cette section nous proposons la résolution numérique d'un problème particulièrement simple : déformation d'une portion de cylindre circulaire sous l'action de pressions internes. Le programme utilisé est général, i.e., opérationnel pour des coques de forme quelconque. Ce programme sera incorporé dans le code MODULEF (cf. PERRONNET [1]).

**Définition de la surface moyenne de la coque considérée** : cf. Figure 3.2-1.

Nous utilisons les notations et les résultats de CIARLET [1]. Pour calculer les $A_{IJ}(x,y)$ il convient d'utiliser les 33 fonctions suivantes (II, 1.3-4) calculées à l'aide des dérivées partielles premières, secondes et troisièmes de la carte $\vec{\phi}(x,y) = R \cos x \vec{e}_1 + R \sin x \vec{e}_2 + y$. Il vient :
<figure><img src="image_1.png" /><figcaption>**Figure 3.2-1**

**Définition de la surface moyenne de la coque**</figcaption></figure>

$$
(3.2-1) \quad \left\{
\begin{array}{l}
a_{11} = R^2, \ a_{12} = a_{21} = 0, \ a_{22} = 1, \ a = R^2, \\
a^{11} = \frac{1}{R^2}, \ a^{12} = a^{21} = 0, \ a^{22} = 1, \\
b_{11} = -R, \ b_{12} = b_{21} = b_{22} = 0, \\
b_1^1 = -\frac{1}{R}, \ b_2^1 = b_1^2 = b_2^2 = 0, \\
b^{11} = -\frac{1}{R^3}, \ b^{12} = b^{21} = b^{22} = 0, \\
c_1^1 = \frac{1}{R^2}, \ c_1^2 = c_2^1 = c_2^2 = 0, \\
\Gamma_{\beta\gamma}^{\alpha} = 0, \ b^{\alpha\beta}|_{\gamma} = 0.
\end{array}
\right.
$$

Calcul explicite de la forme bilinéaire $a(\vec{u}, \vec{v})$ :

Suivant (II, (1.2-1)) la forme bilinéaire $a(\vec{u}, \vec{v})$ s'écrit

$$
(3.2-2) \qquad a(\vec{u}, \vec{v}) = \int_{\Omega} t_{U}[A_{IJ}] V \, dxdy.
$$

Nous allons calculer explicitement la matrice $A_{IJ}$. Suivant (II, (1.3-2))
on obtient :
$$
(3.2-3) \quad \left\{
\begin{array}{l}
\Lambda_1^1 = [0 \quad \frac{1}{R^2} \quad 0 \quad 0 \quad 0 \quad 0 \quad \frac{1}{R} \quad 0 \quad 0 \quad 0 \quad 0 \quad 0], \\
\Lambda_1^2 = [0 \quad 0 \quad \frac{1}{2} \quad 0 \quad \frac{1}{2} \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0], \\
\Lambda_2^1 = [0 \quad 0 \quad \frac{1}{2R^2} \quad 0 \quad \frac{1}{2R^2} \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0], \\
\Lambda_2^2 = [0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 1 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0].
\end{array}
\right.
$$

Suivant (II, (1.3-5)) il vient

$$
(3.2-4) \quad \left\{
\begin{array}{l}
M_1^1 = [0 - \frac{2}{R^3} \quad 0 \quad 0 \quad 0 \quad 0 - \frac{1}{R^2} \quad 0 \quad 0 \quad \frac{1}{R^2} \quad 0 \quad 0], \\
M_1^2 = [0 \quad 0 - \frac{1}{R} \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 1 \quad 0], \\
M_2^1 = [0 \quad 0 - \frac{1}{R^3} \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad \frac{1}{R^2} \quad 0], \\
M_2^2 = [0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 0 \quad 1].
\end{array}
\right.
$$

Alors

$$
(1-\nu)^t \Lambda_{\alpha\beta}^{\beta\alpha} + \nu^t \Lambda_{\alpha\beta}^{\alpha\beta} = \begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \frac{1}{R^4} & 0 & 0 & 0 & \frac{\nu}{R^2} & \frac{1}{R^3} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-\nu}{2R^2} & 0 & \frac{1-\nu}{2R^2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-\nu}{2R^2} & 0 & \frac{1-\nu}{2R^2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & \frac{\nu}{R^2} & 0 & 0 & 0 & 1 & \frac{\nu}{R} & 0 & 0 & 0 & 0 & 0 \\
0 & \frac{1}{R^3} & 0 & 0 & 0 & \frac{\nu}{R} & \frac{1}{R^2} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}
$$
<table><tbody><tr><td>(1-ν) t<sup>β</sup><sub>M<sup>α</sup>αβ</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>+ ν t<sup>α</sup><sub>M<sup>α</sup>αβ</sub> =</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>-1/4</td><td>0</td><td>-ν/2</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>0</td><td>-2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>ν/2</td><td>0</td><td>0</td><td>1</td></tr><tr><td></td><td>0</td><td>0</td><td>-2(1-ν)/R<sup>3</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2(1-ν)/R<sup>2</sup></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>0</td><td>-2ν/R<sup>3</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>-ν/R<sup>2</sup></td><td>0</td><td>0</td><td>ν/R<sup>2</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

Alors avec (II,(1.3-8)) on obtient :
**Formulation variationnelle du problème :**

Suivant (II, (1.1-9), (1.2-1), (1.2-5)) le problème admet la formu-
lation variationnelle suivante :

$$
(3.2-6) \quad \int_{\Omega} t_U [A_{IJ}] V \, dxdy = \int_{\Omega} t_F V \, dxdy, \quad \nabla_{\vec{V}} \in \nabla_{\vec{V}}'' .
$$

Comme la coque est encastrée sur son bord, cette équation est (formel-
lement) équivalente au système d'équations suivant :

$$
(3.2-7) \quad \left\{
\begin{aligned}
& - \frac{\partial}{\partial x} \left[ \frac{Ee}{1-\nu^2} \left\{ \frac{1}{R^4} \left(1+\frac{e^2}{3R^2}\right) u_{1,x} + \frac{\nu}{R^2} u_{2,y} + \frac{1}{R^3} \left(1+\frac{e^2}{6R^2}\right) u_3 - \frac{e^2}{6R^5} u_{3,xx} - \frac{\nu e^2}{6R^3} u_{3,yy} \right\} \right] \\
& - \frac{\partial}{\partial y} \left[ \frac{Ee}{1+\nu} \left\{ \frac{1}{2R^2} \left(1+\frac{e^2}{3R^2}\right) u_{1,y} + \frac{1}{2R^2} u_{2,x} - \frac{e^2}{6R^3} u_{3,xy} \right\} \right] = f^1
\end{aligned}
\right.
$$

$$
(3.2-8) \quad \left\{
\begin{aligned}
& - \frac{\partial}{\partial x} \left[ \frac{Ee}{1+\nu} \frac{1}{2R^2} (u_1, y + u_2, x) \right] - \frac{\partial}{\partial y} \left[ \frac{Ee}{1-\nu^2} \left\{ \frac{\nu}{R^2} u_{1,x} + u_{2,y} + \frac{\nu}{R} u_3 \right\} \right] = f^2 \\
& \frac{Ee}{1-\nu^2} \left[ \frac{1}{R^3} \left(1 + \frac{e^2}{6R^2}\right) u_{1,x} + \frac{\nu}{R} u_{2,y} + \frac{1}{R^2} \left(1 + \frac{e^2}{12R^2}\right) u_3 - \frac{e^2}{12R^4} u_{3,xx} - \frac{\nu e^2}{12R^2} u_{3,yy} \right]
\end{aligned}
\right.
$$

$$
(3.2-9) \quad \left\{
\begin{aligned}
& + \frac{\partial^2}{\partial x^2} \left[ \frac{Ee^3}{12R^5(1-\nu^2)} \left\{ -2u_{1,x} - Ru_3 + Ru_{3,xx} + \nu R^3 u_{3,yy} \right\} \right] \\
& + \frac{\partial^2}{\partial x \partial y} \left[ \frac{Ee^3}{6R^3(1+\nu)} \left\{ -u_{1,y} + Ru_{3,xy} \right\} \right] \\
& + \frac{\partial^2}{\partial y^2} \left[ \frac{Ee^3}{12R^3(1-\nu^2)} \left\{ -2\nu u_{1,x} - \nu Ru_3 + \nu Ru_{3,xx} + R^3 u_{3,yy} \right\} \right] = f^3.
\end{aligned}
\right.
$$

Rappelons que $f^i$, $i=1,2,3$, sont les composantes contravariantes du champ de pression $\vec{f}$ sur la base locale $\vec{a}_i$, i.e.,

$$
\vec{f} = f^i \vec{a}_i.
$$

**Exemple 3.2-1** (E, e, v sont supposés constants, 0≤α<β≤2π)

$$
(3.2-10) \left\{
\begin{aligned}
f^1 &= 0, \\
f^2 &= \frac{4Ee\nu}{(1-\nu^2)R} (H^2-y^2)y, \\
f^3 &= \frac{Ee}{12R^2(1-\nu^2)} \left\{ \left(12 + \frac{e^2}{R^2}\right) (H^2-y^2)^2 + 8ve^2(H^2-3y^2) + 24e^2R^2 \right\}.
\end{aligned}
\right.
$$
La solution exacte est alors

$$ (3.2-11) \begin{cases} u_1 = u_2 = 0, \\ u_3(x,y) = (H^2-y^2)^2 \end{cases} $$

Ce problème correspond au cas où la coque est encastrée sur les bords $y = \pm H$. Les bords $x = \alpha$ et $x = \beta$ sont des bords fictifs sur lesquels on applique des conditions de symétrie liées à la symétrie du champ de pression ($\vec{f}$ est indépendant de $x$ et $f^1 = 0$). La prise en compte de ces conditions aux limites et de ces conditions de symétrie s'effectue comme suit :

### (i) Conditions aux limites (bords $y = \pm H$)

Sur les bords $y = \pm H$, les seuls degrés de liberté non nuls sont obtenus par analogie avec les Figures 2.1-2 et 2.2-2. La prise en compte des conditions aux limites est résumée dans la Figure 3.2-2.

<table><caption>Figure 3.2-2
Prise en compte des conditions aux limites</caption><thead><tr><th rowspan="2">Bords</th><th colspan="6">Sommet d'un triangle</th><th>Milieu d'un côté</th></tr><tr><th>v</th><th>v<sub>,x</sub></th><th>v<sub>,y</sub></th><th>v<sub>,xx</sub></th><th>v<sub>,xy</sub></th><th>v<sub>,yy</sub></th><th>v<sub>,n</sub></th></tr></thead><tbody><tr><td rowspan="3">y = ± H</td><td>v<sub>1</sub></td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>v<sub>2</sub></td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>v<sub>3</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td></tr></tbody></table>

### (ii) Conditions de symétrie (bords $x = \alpha, y = \beta$)

Le champ de pression étant indépendant de $x$, il en va de même du champ de déplacement $\vec{u}$. De plus $f^1 = 0$ entraîne $u_1 = 0$ en tous les points de la coque (rappelons que $\vec{u} = u_i \vec{a}^i$), d'où $\vec{u}(x,y) = u_2(y)\vec{a}^2 + u_3(y)\vec{a}^3$. Alors la prise en compte des conditions de symétrie comporte les deux parties suivantes

* **Suppressions de certains degrés de liberté, conformément à la Figure 3.2-3 ;**

** **Formation de la matrice de rigidité globale en assemblant les matrices de rigidité élémentaires des triangles recouvrant $\Omega$ mais également des triangles extérieurs voisins (en pointillés sur la Figure 3.2-4).**
<table><caption>Figure 3.2-3
Suppressions liées aux conditions de symétrie</caption><thead><tr><th rowspan="2">Bords</th><th colspan="6">Sommet d'un triangle</th><th>Milieu d'un côté</th></tr><tr><th>v</th><th>v,<sub>x</sub></th><th>v,<sub>y</sub></th><th>v,<sub>xx</sub></th><th>v,<sub>xy</sub></th><th>v,<sub>yy</sub></th><th>v,<sub>n</sub></th></tr></thead><tbody><tr><td>x = α</td><td>v1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>x = β</td><td>v2</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>v3</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 3.2-4
Pour la formation de la matrice de rigidité globale il faut prendre en compte les triangles situés dans Ω et les triangles voisins (en pointillés).</figcaption></figure>

Remarque 3.2-1 : Dans l'exemple 3.2-1 le champ de déplacement $\vec{u}$ est indépendant de $x$. Pour le déterminer on peut se limiter à un domaine $\Omega$ pour lequel $\beta = \alpha+h$ de telle sorte que le maillage se limite à un seul triangle dans la direction des $x$ (cf. Figure 3.2-5).
<figure><img src="image_3.png" /><figcaption>Figure 3.2-5</figcaption></figure>

Lorsque le champ de pressions est indépendant de x les considéra-
tions de symétrie permettent de se limiter à la "bande" [α, α+h].

**Exemple 3.2-2** (E,e,v sont supposés constants, 0≤α<β≤2π)

$$
(3.2-12) \quad \begin{cases} f^1 = 0, \\ f^2 = \frac{2Ee}{1-v^2}, \\ f^3 = - \frac{2Eev}{R(1-v^2)} y. \end{cases}
$$

La solution exacte est alors

$$
(3.2-13) \left\{
\begin{array}{l}
u_1 = u_3 = 0, \\
u_2 = H^2 - y^2.
\end{array}
\right.
$$

Pour la prise en compte des conditions aux limites et des conditions
de symétrie, nous renvoyons à l'Exemple 3.2-1 et à la Remarque 3.2-1.

**Exemple 3.2-3** (E,e,v sont supposés constants, 0≤α<β≤2π)

$$
(3.2-14) \left\{
\begin{array}{l}
f^1 = \frac{Ee}{(1+\nu)R^2} \left(1 + \frac{e^2}{3R^2}\right) \\
f^2 = f^3 = 0
\end{array}
\right.
$$
La solution exacte est alors

$$ (3.2-15) \begin{cases} u_1 = H^2-y^2 \\ u_2 = u_3 = 0 \end{cases} $$

Dans ce cas, la prise en compte des conditions aux limites s'effectue comme indiqué sur la Figure 3.2-2. Pour la prise en compte des conditions de symétrie, les Figures 3.2-4 et 3.2-5 restent valables et il convient de remplacer la Figure 3.2-3 par la suivante (u̇ est indépendant de x) :

<table><caption>Figure 3.2-6
Suppressions liées aux conditions de symétrie dans l'Exemple 3.2-3.</caption><thead><tr><th rowspan="2">Bords</th><th colspan="6">Sommet d'un triangle</th><th rowspan="2">Milieu d'un côté</th></tr><tr><th>v</th><th>v,<sub>x</sub></th><th>v,<sub>y</sub></th><th>v,<sub>xx</sub></th><th>v,<sub>xy</sub></th><th>v,<sub>yy</sub></th></tr></thead><tbody><tr><td>x = α</td><td>v<sub>1</sub></td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>x = β</td><td>v<sub>2</sub></td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td></td><td>v<sub>3</sub></td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td></tr></tbody></table>

## Résultats numériques.

Suivant la Section 3.1, nous avons implémenté la méthode d'approximation par éléments finis de type ARGYRIS pour la résolution de problèmes de coques généraux, puis nous avons testé ce programme sur les Exemples 3.2-1 à 3.2-3. Les résultats s'avèrent excellents, ce qui était prévisible si l'on observe que les solutions exactes appartiennent à l'espace des fonctions interpolantes.
BIBLIOGRAPHIE

L'essentiel de la bibliographie figure dans II et dans BERNADOU-BOISSERIE [1].
Nous reprenons ici les références citées dans ce rapport.

BERNADOU M.,

[1] Méthodes conformes d'éléments finis avec intégration numérique
pour des problèmes de coques. Rapport IRIA-LABORIA N° 245, Août 1977.

[2] Sur l'analyse numérique du modèle linéaire de coques minces de
W.T. KOITER. Thèse d'Etat, Paris, Juin 1978.

BERNADOU M., BOISSERIE J.M.,

[1] Implémentation de l'élément fini d'ARGYRIS. Exemples. Rapport
IRIA-LABORIA N° 301, Mai 1978.

CIARLET P.G.,

[1] La méthode des éléments finis appliquée aux coques.
Rapport IRIA-LABORIA 113 (1975).

KOITER W.T.,

[1] On the nonlinear theory of thin elastic shells. Proc. Kon. Ned.
Akad. Wetensch. B69 (1966), 1-54.

[2] On the foundations of the linear theory of thin elastic shells.
Proc. Kon. Ned. Akad. Wetensch., B73, (1970), 169-195.

PERRONNET A.,

[1] The Club MODULEF, a library of subroutines for finite element analysis. Proceedings of the "Third International Symposium on Computing Methods in Applied Sciences and Engineering", Versailles (1977).
ANNEXE A

On trouvera ci-après la liste du programme principal relatif à la résolution
de problèmes de coques.
IMPLICIT REAL *8(A-H),REAL *8(O-Z)
COMMON/BURA/IRA,IRA1,IRA2,IRB,IRB1,IRB2,IQUVER,NPOIN,IRAL
COMMON/BURB/NINCO,NINCOA,MINCO,MJNCO
COMMON/BURR/H1,H2,R,XZER,YZER
COMMON/BERNAD/YOUNG,POISS,EPAISS
COMMON/LECBER/ILEC,IVEBER
DIMENSION IER(441),IFR(441),IKR(1046),IHR(1046)
DIMENSION A(15000)
DATA DEUPI/6.2831853071796/
DATA UN/1.DO/
IDIMA=15000
YOUNG=100.DO
XZER=0.DO
YZER=-0.5DO
POISS=0.3DO
EPAISS=0.01
R=1.DO
ILEC=0
ILEC=1
IQUVER=1
IQUVER=0
IRA=7
IRB=5
LAS=1 POTENTIEL
LAS=2 PLAQUE
LAS=3 COQUE COMPLETE
LAS=3
IVEBER=LAS-1
NINCO=1+2*(IVEBER/2)
MINCO=NCOCO(NINCO,MJNCO)
NINCOA SERT A COMPACTER DE L'INFORMATION
NINCOA=100**MINCO
IF(IQUVER.EQ.1) DEUPI=UN
IF(IQUVER.EQ.1)YZER=XZER
H1=DEUPI/DFLOAT(IRA-1)
H2=UN/DFLOAT(IRB-1)
IF(IVEBER.EQ.2) CALL INGDIF
CALL P5BER1(A,IER,IFR,IHR,IKR,IDIMA,IRA*IRB)
1 CONTINUE
STOP
END
ANNIMP. B

Dans cette annexe nous rassemblons les sous-programmes relatifs aux données
des exemples traités dans la section 3.2, I.a.,

;

FFF

TEMP

XDO-XD1-XD2-XD3
FUNCTION FFF(X,Y,INCO)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

C TEST DU CYLINDRE (THESE DE BERNADOU) FORMULES (3-2-10)

COMMON/BERNAD/ YOUNG,POISS,EPAISS

COMMON/BURR/H1,H2,R

DATA ZERO,UN/0.00,1.00/

DATA TROIS,QUATRE,HUIT,DOUZE/3.00,4.00,8.00,12.00/

DATA H/0.500/

UU(G,Z)=G+G-Z*Z

FFF=ZERO

POI2=UN-POISS**2

IF(INCO-2)=1,2,3

1 RETURN

2 FFF=QUATRE*YOUNG*EPAISS*(POISS/POI2)*UU(H,Y)*(Y/R)

RETURN

3 FFF={YOUNG*EPAISS/(DOUZE*R*R*POI2)}*

1 (((DOUZE+{EPAISS/R})**2)*UU(H,Y)**2)

2 +HUIT*POISS*EPAISS**2*(UU(H,Y)-2.00*Y*Y)

3 +TROIS*HUIT*(EPAISS*R)**2)

RETURN

END

FUNCTION TEMP(XA,INC,INCO)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

C TEST DU CYLINDRE (THESE DE BERNADOU) FORMULES (3-2-11)

COMMON/LECBER/ILEC,IVEBER

DIMENSION XA(2)

DATA H,DEUX,QUATRE/0.500,2.00,4.00/

DATA ZERO/0.00/

UU(G,Z)=G+G-Z*Z

TEMP=ZERO

IF(INCO.LT.3) RETURN

K1=K1TEMP(INC)

GO TO(1,2,3,4,5,6),K1

1 TEMP=UU(H,XA(2))**2

2 RETURN

3 TEMP=-QUATRE*XA(2)*UU(H,XA(2))

4 CONTINUE

5 RETURN

6 TEMP=-QUATRE*(UU(H,XA(2))-DEUX*XA(2)**2)

RETURN

END
FUNCTION XDO(A,B,IVER)

IMPLICIT REAL*8(A-H),REAL*8(O-Z)

GENERIC

COMMON/BURR/H1,H2,R,XZER,YZER

DATA ZERO,UN/O,DO,1,DO/

IF(IVER-2) 10,20,30

10 XDO=R*COS(A)

RETURN

20 XDO=R*SIN(A)

RETURN
30 XDO=B

RETURN

ENTRY XD1(A,B,IVER)

IX=MOD(IVER,10)

IDX=IVER/10

1 IF(IX-2) 11,21,31

11 XD1=-R*SIN(A)

RETURN

21 XD1= R*COS(A)

RETURN
31 XD1=ZERO

RETURN
2 IF(IX-2) 12,22,32

12 XD1 = ZERO

RETURN
22 XD1 = ZERO

RETURN
32 XD1 = UN

RETURN

ENTRY XD2(A,B,IVER)

IX=MOD(IVER,10)

IDX=IVER/10

3 IF(IX-2) 13,23,33

13 XD2=- R*COS(A)

RETURN

23 XD2=- R*SIN(A)

RETURN

33 XD2=ZERO

RETURN

4 IF(IX-2) 14,24,34
14

RETURN

24

RETURN

34

RETURN

5

IF(IX-2) 15,25,35

15

RETURN

25

RETURN

35

RETURN

ENTRY XD3(A,B,IVER)

IX=MOD(IVER,10)

IDX=IVER/10

GO TO (6,7,8,9),IDX

6

IF(IX-2) 16,26,36

16

RETURN

26

RETURN

36

RETURN

7

IF(IX-2) 17,27,37

17

RETURN

27

RETURN

37

RETURN

8

IF(IX-2) 18,28,38

18

RETURN

28

RETURN

38

RETURN

9

IF(IX-2) 19,29,39

19

RETURN

29

RETURN

39

RETURN

END

XD2=ZERO

XD2=7ERO

XD2=ZERO

XD2=ZERO

XD2=ZERO

XD2=ZE#3

XD3=R*SIN(A)

XD3=-R*COS(A)

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZERO

XD3=ZER#

XD3=ZER#
# ANNEXE C

Dans cette annexe nous rassemblons les sous-programmes relatifs au problème de *coque*. Nous renvoyons à BERNADOU-BOISSERIE [1] pour la liste des sous-programmes relatifs à la triangulation, à l'interpolation, à l'intégration, à certaines parties de l'assemblage, à la résolution ainsi que pour les sous-programmes de service.

## C-1 - MODULE GEOMETRIE DIFFERENTIELLE.

Ce module comporte :

- FEGEO
- FEPHI
- INGDIF
- MATA-MATB
- MAT2
- METGOS
- MIXTE
- PROSCA
- PROVE1-PROVE2
- SORTIE
- STOFFE
## C-2 - MODULE RIGIDITE ELEMENTAIRE.

Ce module comporte :

ESCLA1

ESCLA2

FEAIJ

FEAIJ1

MACOQ1

MACOQ2

MACOQ3

## C-3 - MODULE SECOND MEMBRE ELEMENTAIRE.

Ce module comporte :

K1TEMP

MBCOQ3
FUNCTION ESCLA1(Phi,ACONTR,BMIXT,BCONTR,TLFE,RACINA)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

COMMON/BURF/IALF,IBET,IGAM

DIMENSION Phi(3,2),ACONTR(3),BCONTR(3),BMIXT(4),TOFFE(6)

DIMENSION IA(2,2),IB(2,2)

DATA IA/1,2,2,3/

DATA IB/1,2,3,4/

DATA ZERO/0,00/

ROM=ZERO

SOM=ZERO

TOM=ZERO

DO 1 ILAM=1,2

DO 1 INUU=1,2

I1=IA(IALF,ILAM)

I2=IA(1BET,1NUU)

I3=IGAM+IA(ILAM,INUU)+8

I4=IA(ILAM,INUU)+3

I5=IGAM+7

A1=ACONTR(I1)*ACONTR(I2)

A2=PROSCA(3,Phi(1,7),Phi(1,I3))

A3=PROSCA(3,Phi(1,I5),Phi(1,I4))

SOM=SOM+A1*(A2+A3)

1 CONTINUE

DO 2 ILAM=1,2

DO 2 INUU=1,2

I1=IA(IALF,ILAM)

I2=IA(1BET,1LAM)

I3=IB(1BET,1NUU)

I4=I3(IALF,INUU)

J1=IA(ILAM,IGAM)

I5=2*J1+(INUU-2)

A1=ACONTR(I1)*BMIXT(I3)

A2=ACONTR(I2)*BMIXT(I4)

TOM=TOM+TOFFE(I5)*(A1+A2)

2 CONTINUE

IO=IA(IALF,IBET)

DO 3 INUU=1,2

J1=IA(INUU,IGAM)

I1=2*J1+(INUU-2)

ROM=ROM+BCONTR(IO)*TOFFE(I1)

3 CONTINUE

ESCLA1=SOM*RACINA-TOM-ROM

RETURN

END

SUBROUTINE ESCLA2(ACONTR,BCONTR,BMIXT,CAB,COUR,TOF,AL,AM)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

COMMON/BURF/IALF,IBET,IGAM

DIMENSION ACONTR(3),BCONTR(3),BMIXT(4),CAB(4),TOF(2,3),C IR(2,3)

DIMENSION AL(4,12),AM(4,12),DELT(4),IA(2,2)

DATA DELT/1,0,0,0,0,1,00/

DATA IA/1,2,2,3/

DATA ZERO,DEMI,UN,DEUX/0,0,0,500,1,000,0,000

DEFINITION DE LA......LIEE

ILIG=IALF+2*(1BET-1)

DO 1 J=1,12

AL(ILIG,J)=ZER

1 AM(ILIG,J)=ZER

I1=ILIG

AL(ILIG,7)=-BM1 T1)

AM(ILIG,7)=-CAB1
DO 2 K=1,2

J=4*K-2

J2=8+2*K

I1=IA(K,IALF)

I2=IBET+2*(K-1)

AL(ILIG,J)=ACONTR(I1)*DELT(I2)

J0=2*K+1

DO 3 ILAM=1,2

I1=IA(IALF,ILAM)

I2=ILAM+2*(IBET-1)

3 AL(ILIG,J0)=AL(ILIG,J0)+DEMI*(UN-DELT(I2))*ACONTR(I1)

J1=3*K-2

DO 4 INUU=1,2

I1=IA(IALF,INUU)

I2=IA(IBET,INUU)

4 AL(ILIG,J1)=AL(ILIG,J1)-ACONTR(I1)*TOF(K,I2)

2 CONTINUE

C

CALCUL DES......... AM

AM(ILIG,10)=AL(ILIG,2)

AM(ILIG,11)=DEUX*AL(ILIG,3)

AM(ILIG,12)=AL(ILIG,6)

DO 20 K=1,2

J=1+3*(K-1)

J0=7+K

J1=8+2*K

J2=3*K-2

I1=IA(K,IALF)

AM(ILIG,J)=COUR(1BET,I1)

AM(ILIG,J0)=AL(ILIG,J2)

DO 21 ILAM=1,2

I2=IA(IALF,ILAM)

I3=IA(ILAM,IBET)

I4=ILAM+2*(IBET-1)

BMIXTA=BMIXT(I4)

DO 22 INUU=1,2

I20=IA(IALF,INUU)

I22=IA(ILAM,INUU)

AM(ILIG,J)=AM(ILIG,J)-ACONTR(I20)*BMIXTA*TOF(K,I22)

22 CONTINUE

AM(ILIG,J)=AM(ILIG,J)-BCONTR(I2)*TOF(K,I3)

21 CONTINUE

20 CONTINUE

DO 30 K=1,2

I3=IBET+2*(K-1)

I4=IA(K,IALF)

DO 30 L=1,2

J=K+3*L-2

I1=IA(IALF,L)

I2=L+2*(IBET-1)

AM(ILIG,J)=BCONTR(I1)*DELT(I3)+A(14)

30 CONTINUE

RETURN

END
SUBROUTINE FEAIJ
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
COMMON/BURM/IADS,IACOV,IACOT,IBCOV,IBCOT,IBMIX,ITOFF,ICOUR,ICAB,
IAL,IAM
COMMON/BURN/                      ALPHA,BETA
COMMON/BURT/IAIJ1
COMMON/BURU/AMESUR
COMMON/LECBER/ILEC,IVBER
DIMENSION A(129),PHI(3,13)

CALL NULMAT(IAG1,1,A)
CALL FEPHI(ALPHA,BETA,PHI)
CALL FEGEO(              PHI,
1           A(IACOV),A(IACOT),A(IADS),
2           A(IBCOV),A(IBCOT),A(IBMIX),
3           A(ITOFF),A(ICAB),A(ICOUR),A(IAL),A(IAM))
AMESUR=A(IADS)
IF(IAIJ1.EQ.2) RETURN
IF(ILEC.NE.50) GO TO 1
CALL SORTIE(          PHI,
1           A(IACOV),A(IACOT),A(IADS),
2           A(IBCOV),A(IBCOT),A(IBMIX),
3           A(ITOFF),A(ICAB),A(ICOUR),A(IAL),A(IAM))
1 CONTINUE
CALL FEAIJ1(A(IAL),A(IAM))

END

SUBROUTINE FEAIJ1(AL,AM)
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
LOGICAL VRAI
GENERIC
COMMON/BERNAC/APOIDS
COMMON/BERNAD/YOUNG,POISS,EPAISS
COMMON/BERNAE/LOGG(12,12),AIJ(78)
COMMON/BURS/IELEM,INT
COMMON/BURU/AMESUR
COMMON/LECBER/ILEC,IVBER
DIMENSION G(12,12),AL(4,12),AM(4,12),COE(4)
DATA ZERO,ZERO1/0.00,1.0-16/
DATA UN,DOUZE/1.00,12.00/
1100 FORMAT(1X,11(G10.3,1X),G10.3)
1101 FORMAT(1H1,26HEDITION DE LA MATRICE AIJ/1X,26(1H*)/1H-)
1102 FORMAT(1X,11(I3,8X),I3)
C
FONCTION ARITHMETIQUE DEFINIE AU DEBUT DU S/P
IL(IA,IB)=IA+2*(IB-1)
VRAI=.FALSE.
IF((ILEC.NE.5).OR.(INT.NE.1).OR.(IELEM.GT.02)) VRAI=.TRUE.
COEFF=EPAISS*EPAISS/DOUZE
COEG=YOUNG*EPAISS*APOIDS*SQRT(AMESUR)/(UN-POISS*POISS)
COE(1)=UN-POISS
COE(2)=POISS
COE(3)=COEFF+(UN-POISS)
COE(4)=COEFF+POISS
DO 100 I=1,4

100 COE(I)=COEFG*COE(I)
CALL NULMAT(12,12,G)
DO 1 IALPHA=1,2
DO 1 IBETA=1,2
IL1=IL(IBETA,IALPHA)
IL2=IL(IALPHA,IBETA)
IL3=IL(IALPHA,IALPHA)
IL4=IL(IBETA,IBETA)
DO 2 I=1,12
I1=I

C SEULE LA MATRICE AIJ SUPERIEURE EST CONSTITUEE
DO 2 J=I1,12
LOGG(I,J)=0

1 +COE(1) *AL(IL1,I)*AL(IL2,J)

2 +COE(2) *AL(IL3,I)*AL(IL4,J)

3 +COE(3) *AM(IL1,I)*AM(IL2,J)

4 +COE(4) *AM(IL3,I)*AM(IL4,J)

2 CONTINUE

1 CONTINUE
IF(.NOT.VRAI) PRINT 1101
ICOMPT=0
GMAX=ABS(G(1,1))
DO 3 I=1,12
I1=I

DO 5 J=I1,12
IF(ABS(G(J,I)).GT.GMAX) GMAX=ABS(G(J,I))

5 G(J,I)=G(I,J)

4 CONTINUE
DO 6 J=1,12
GTEST=ABS(G(I,J))
IF(GMAX.NE.ZERO) GTEST=GTEST/GMAX
IF(GTEST.GT.ZERO1) GO TO 8

C OPERATION DE NETTOYAGE DE LA MATRICE G
G(I,J)=ZERO
GO TO 6

8 CONTINUE
IF(I.GT.J) GO TO 6
ICOMPT=ICOMPT+1
LOGG(I,J)=ICOMPT
LOGG(J,I)=ICOMPT
AIJ(ICOMPT)=G(I,J)

6 CONTINUE
IF(VRAI) GO TO 3
PRINT 1100,(G(I,K),K=1,12)

3 CONTINUE
IF(.NOT.VRAI) CALL LIVEC(AIJ,ICOMPT,ICOMPT+2000)
IF(VRAI) RETURN
CALL AVANPA(-1)
DO 9 I=1,12
PRINT 1102,(LOGG(I,K),K=1,12)

9 CONTINUE
RETURN
END
SUBROUTINE FEGEO(PHI,ACOVAR,ACONTR,ADS,BCOVAR,BCONTR,BMIXT,TOF,CAB
1,COUR,AL,AM)

IMPLICIT REAL*8(A-H),REAL*8(O-Z)

GENERIC

COMMON/BURF/IALF,IBET,IGAM

DIMENSION PHI(3,2),ACOVAR(3),ACONTR(3),BCOVAR(3),BCONTR(3),IA(3,2)

DIMENSION TOF(2,3),BMIXT(4),COUR(2,3),CAB(2,2),IB(2,2)

1 ,AL(4,2),AM(4,2)

DATA IA/2,2,3,2,3,3/

DATA IB/1,0,2,3/

DATA UN/1.00/

DO 1 I=1,3

I1=IA(I,1)

I2=IA(I,2)

ACOVAR(I)=PROSCA(3,PHI(1,I1),PHI(1,I2))

1 CONTINUE

CALL MAT2(ACOVAR,ACONTR,ADS)

RACINA=UN/SQRT(ADS)

DO 2 I=1,3

I1=I+3

BCOVAR(I)=PROSCA(3,PHI(1,7),PHI(1,I1))*RACINA

2 CONTINUE

CALL MIXTE(ACONTR,BCOVAR,BMIXT,BCONTR,CAB)

CALL STOFFE (PHI,ACONTR,TOF)

DO 5 IGAM=1,2

DO 5 IALF=1,2

DO 5 IBET=1,2

I1=IB(IALF,IBET)

IF(I1.EQ.0) GO TO 5

COUR(IGAM,I1)=ESCLA1(PHI,ACONTR,BMIXT,BCONTR,TOF,RACINA)

5 CONTINUE

DO 10 IALF=1,2

DO 10 IBET=1,2

CALL ESCLA2(ACONTR,BCONTR,BMIXT,CAB,COUR,TOF,AL,AM)

10 CONTINUE

RETURN

END
SUBROUTINE FEPHI(A,B,PHI)
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
DIMENSION PHI(3,2)
DO 1 I=1,3
I1=I
PHI(I,1)=XDO(A,B,I1)
DO 2 J=2,3
I2=I+10*(J-1)
2 PHI(I,J)=XD1(A,B,I2)
DO 3 J=4,6
I3=I+10*(J-3)
3 PHI(I,J)=XD2(A,B,I3)
DO 4 J=10,13
I4=I+10*(J-9)
4 PHI(I,J)=XD3(A,B,I4)
1 CONTINUE
I2=3
DO 5 I=1,3
I0=(I/3)
I1=2*I-I0
I3=I+6
CALL PROVE1(PHI(1,I1),PHI(1,I2),PHI(1,I3))
5 CONTINUE
I1=2
DO 7 I=1,2

I2=I+4
I3=7+I
CALL PROVE2(PHI(1,I1),PHI(1,I2),PHI(1,I3))
7 CONTINUE
RETURN
END

SUBROUTINE INGDIF
COMMON/BURH/IADS,IACOV,IACOT,IBCOV,IBCOT,IBMIX,ITOFF,ICOUR,ICAB,
IAL,IAM,IAG1
IADS=1
IACOV=2
IACOT=IACOV+3
IBCOV=IACOT+3
IBCOT=IBCOV+3
IBMIX=IBCOT+3
ITOFF=IBMIX+4
ICAB=ITOFF+6
ICOUR=ICAB+4
IAL=ICOUR+6
IAM=IAL+48
IAG1=IAM+47
RETURN
END

FUNCTION K1TEMP(INC)
COMMON/BURH/NTYP1,NTYP2,ITYPZ(6),ITYP2(6)
IVER=INC/10000
IF(IVER.EQ.0) K1TEMP=ITYP2(INC)
IF(IVER.EQ.1) K1TEMP=MOD(INC,10000)
RETURN
END
SUBROUTINE MACOQ1(AUX,C,ITRANS,JTRANS,NAI,NAJ,NC)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

LOGICAL VRAI

AUXILIAIRE DE MACOQ3 ... EXPLOITE AIJ

DIMENSION AUX(21,6),C(3,21),C(NC,2)

COMMON/BERNAE/LOGG(12,12),AIJ(78)

COMMON/LEGBER/ILEC,IVEBER

AIJ REPRESENTE LE NOYAU DE L OPERATEUR

AUX EST UNE DONNEE DE LA MACRO-INSTRUCTION QUI VA PRODUIRE C
ITRANS,JTRANS SERVENT A DEFINIR LA PARTIE DE AIJ

*(*)*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*task:ocr>
SUBROUTINE MACOO3(AMAT,X)

IMPLICIT REAL *8(A-H),REAL *8(O-Y)

GENERIC

C
ELEMENT DE COQUE A EPAISSEUR CONSTANTE DANS R3
COMMON/BERNAA/DELT,DELT2,BG(2,3)
COMMON/BERNAB/ALSOM(3),AN(3,3),AUX(21,6),DA(21,21),GL(3,6),
1
    PSTT(21,10)
COMMON/BERNAC/APOIDS
COMMON/BERNAD/ YOUNG,POISS,EPAISS
COMMON/BURN/XO(2)
COMMON/BURS/IELEM,INT
COMMON/BURT/IAIJ1
COMMON/LECBER/ILEC,IVEBER
DIMENSION AMAT(63,63),X(2,3)

C

C
CALCUL DE LA MATRICE DE RIGIDITE DE COQUE

C

DATA DEMI,UN,DOUZE/0.5DO,1.DO,12.DO/

C
CECI VIA LE COMMUN BURT POUR GERER FEAIJ
I A I J 1=1
NA=63
IF(ILEC.EQ.3) CALL LIVEC(X,6,6)
CALL NULMAT(NA,NA,AMAT)
CALL TRIAGO(X,AN,DA)

C
DELTA=DEUX*SURFACE DE L ELEMENT EN VALEUR ALGEBRIQUE

C
BOUCLE SUR LES NOEUDS D INTEGRATION
DO 1 INT=1,16
CALL INT16(ALSOM,APOIDS,INT)
CALL LAMBX(X,ALSOM,XO)

C
CALCUL DES POIDS EFFECTIFS DU SCHEMA D INTEGRATION
A P O I D S =DEMI*APOIDS*ABS(DELT)

C
CALCUL DES PUISSANCES DE LAMBDA AU POINT D INTEGRATION
CALL PLO(GL,ALSOM)
CALL FEAIJ
DO 30 KVER=1,6
DO 19 L=1,21

19 PSTT(L,KVER)=PL2(L,KVER,GL)
IF(KVER.GT.3) GO TO 30
K1=KVER+6
DO 20 L=1,21

20 PSTT(L,K1)=PL1(L,KVER,GL)
IF(KVER.GT.1) GO TO 30
DO 21 L=1,21

21 PSTT(L,10)=PL(L,0,GL)

30 CONTINUE
DO 3 IVER=1,3
SIGNE=UN

SIGNE=-UN

IF(IVER.EQ.2)

IF(IVER.GT.1) GO TO 6

C
INTERPOLEE DE LA FONCTION
DO 8 L=1,21
8 AUX(L,IVER) =PSTT(L,10)
GO TO 2

6 CONTINUE
J=4-IVER

C INTERPOLEE DE LA DERIVEE PREMIERE

DO 9 L=1,21
AUX(L,IVER) = (BG(J,1)*PSTT(L,7)
1      +BG(J,2)*PSTT(L,8)
2      +BG(J,3)*PSTT(L,9))

3

*SIGNE/ D E L T

9 CONTINUE

2 CONTINUE

M=2-(IVER/3)

N=2-(IVER/2)

SIGN2=SIGNE*DELT2

C INTERPOLEE DE LA DERIVEE SECONDE

DO 4 L=1,21
AUX(L,IVER+3) = (BG(M,1)*BG(N,1) * PSTT(L,1)
1      +BG(M,2)*BG(N,2) * PSTT(L,4)
2      +BG(M,3)*BG(N,3) * PSTT(L,6)
3      +(BG(M,1)*BG(N,2)+BG(M,2)*BG(N,1))*PSTT(L,2)
4      +(BG(M,2)*BG(N,3)+BG(M,3)*BG(N,2))*PSTT(L,5)
5      +(BG(M,3)*BG(N,1)+BG(M,1)*BG(N,3))*PSTT(L,3))
1      *SIGN2

4 CONTINUE

3 CONTINUE

C REALISATION DE (3-1-7)

CALL MACOQ1(AUX,AMAT(01,01),0,0,3,3,NA)

CALL MACOQ1(AUX,AMAT(01,22),0,3,3,3,NA)

CALL MACOQ1(AUX,AMAT(01,43),0,6,3,6,NA)

CALL MACOQ1(AUX,AMAT(22,22),3,3,3,3,NA)

CALL MACOQ1(AUX,AMAT(22,43),3,6,3,6,NA)

CALL MACOQ1(AUX,AMAT(43,43),6,6,6,6,NA)

1 CONTINUE

C (DA) EST GENEREE ET STOCKEE SOUS FORME TRANSPOSEE

DO 40 I=1,3

I1=1+21*(I-1)

I2=I

DO 40 J=I2,3

J1=1+21*(J-1)

CALL MATRIX(20,21,21,21,AMAT(I1,J1),NA,DA,21,AMAT(I1,J1),NA)

40 CONTINUE

C TRANSPOSITION DE DA

CALL MATRIX(0,21,21,21,DA,21,DA,21,DA,21)

DO 41 I=1,3

I1=1+21*(I-1)

I2=I

DO 41 J=I2,3

J1=1+21*(J-1)

CALL MATRIX(20,21,21,21,DA,21,AMAT(I1,J1),NA,AMAT(I1,J1),NA)

41 CONTINUE

C ET MAINTENANT ON COMPLETE AMAT PAR SYMETRIE

C CECI POURRAIT ETRE AMMELORE EN TRAVAILLANT ASSBER
CALL MACOQ2(AMAT,NA)

C EDITION DE LA MATRICE DE RIGIDITE ELEMENTAIRE SI ILEC=3

IF(ILEC.EQ.3) CALL REASY2(AMAT,AMAT,NA,1,NA,1,NA,1)

RETURN

END
SUBROUTINE MATA(A,B,C)
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
DIMENSION A(2,2),B(2,2),C(2,2)
DATA ZERO/O.DO/
DO 1 I=1,2
DO 1 J=1,2
C(I,J)=ZERO
DO 1 K=1,2
C(I,J)=C(I,J)+A(I,K)*B(J,K)
1 CONTINUE
RETURN
ENTRY MATB(A,B,C)
DO 2 I=1,2
DO 2 J=1,2
C(I,J)=ZERO
DO 2 K=1,2
C(I,J)=C(I,J)+A(I,K)*B(K,J)
2 CONTINUE
RETURN
END

SUBROUTINE MAT2(A,B,DET)
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
1100 FORMAT(1X,11HERREUR MAT2/1X,11(1H*))
DIMENSION A(3),B(3)
DATA UN,ZERO/1.DO,O.DO/
DET=A(1)*A(3)-A(2)*A(2)
IF(DET.EQ.ZERO) GO TO 1
DET1=UN/DET
B(1)=A(3)*DET1
B(3)=A(1)*DET1
B(2)=-A(2)*DET1
RETURN
1 CONTINUE
PRINT 1100
CALL LIVEC(A,3,1003)
STOP
END
SUBROUTINE MBCOQ3(B,X)

IMPLICIT REAL *8(A-H),REAL *8(O-Z)

C CALCUL DU SECOND MEMBRE DES PROBLEMES DE COQUES

GENERIC

COMMON/BERNAA/DELT,DELT2,BG(2,3)

COMMON/BERNAB/ALSON(3),AN(3,3),AUX(21,6),DA(21,21),GL(3,6),

1 PSTT(21,10)

COMMON/BURN/XO(2)-

COMMON/BURS/IELEM,INT

COMMON/BURT/IAIJ1

COMMON/BURU/AMESUR

DIMENSION B(63),X(2,3)

DATA DEMI/0.500/

C CECI VIA LE COMMUN BURT POUR GERER FEAIJ

I A I J 1=2

CALL TRIAGO(X,AN,DA)

CALL NULMAT(63,1,B)

DO 1 I=1,16

CALL INT16(ALSON,APOIDS,I)

CALL PLO(GL,ALSON)

CALL LAMBX(X,ALSON,XO)

CALL FEAIJ

C DELTA=DEUX*SURFACE DE L ELEMENT EN VALEUR ALGEBRIQUE

A P O I D S =DEMI*APOIDS*ABS(DELT)

*SQRT(AMESUR)

1 DO 2 L=1,21

2 AUX(L,1)=PL(L,0,GL)

DO 3 INCO=1,3

I1=21*(INCO-1)

DO 3 K=1,21

3 B(I1+K)=B(I1+K)+APOIDS*AUX(K,1)*FFF(XO(1),XO(2),INCO)

1 CONTINUE

DO 4 I=1,3

I1=1+21*(I-1)

CALL MATRIX(20,1,21,21,B(I1),1,DA,21,B(I1),1)

4 CONTINUE

RETURN

END
SUBROUTINE METGOS(PHI,TOF,BCOVAR,ACOV,COUR)
    IMPLICIT REAL*8(A-H),REAL*8(O-Z)
    DIMENSION PHI(3,12),VAR(3),TOCOVA(2,3),TOCOVD(2,3,2),TOF(2,3)
    DIMENSION IA(2,2),IB(2,2,2),EPSI(2,2),COUR(6),ACOV(3)
1100 FORMAT(9X,4HIALF,1X,I1,1X,4HIBET,1X,I1,1X,4HILAM,1X,I1,1X,4HIMUU,
    11X,I1,1X,2HA1,1X,D15.8,1X,2HA2,1X,D15.8,1X,2HA3,1X,D15.8,1X,2HA4,1
    2X,D15.8,1X,3HSOM,1X,D15.8/9X,4(1H*),3X,4(1H*),3X,4(1H*),3X,4(1H*),
    13X,2(1H*),17X,2(1H*),17X,2(1H*),17X,2(1H*),17X,3(1H*))
1101 FORMAT(1X,16HCONTROLE-CODAZZI,20X,2(10X,D15.8)/1X,16(1H*))
    DATA EPSI/0.00,-1.00,1.00,0.00/
    DATA IA/1,2,2,3/
    DATA IB/1,2,2,3,2,3,3,4/
    DATA ZERO/0.00/
    DO 1 IBET=1,2
    DO 2 IMUU=1,2
    I1=IA(1BET,IMUU)
    I3=I1+3
    DO 3 IALF=1,2
    IALF1=IALF+1
    TOCOVA(IALF,I1)=PROSCA(3,PHI(1,IALF1),PHI(1,I3))
    DO 10 ILAM=1,2
    I2=IA(IALF,ILAM)+3
    I4=IB(1BET,IMUU,ILAM)+9
    A1=PROSCA(3,PHI(1,I2),PHI(1,I3))
    A2=PROSCA(3,PHI(1,IALF1),PHI(1,I4))
    TOCOVD(IALF,I1,ILAM)=A1+A2
10 CONTINUE
    3 CONTINUE
    2 CONTINUE
    1 CONTINUE
    DO 4 IALF=1,2
    DO 4 IBET=1,2
    DO 4 ILAM=1,2
    DO 4 IMUU=1,2
    A1=EPSI(IALF,IBET)*EPSI(ILAM,IMUU)
    I1=IA(IALF,IMUU)
    I2=IA(1BET,ILAM)
    I3=IA(1BET,IMUU)
    A2=BCOVAR(I1)*BCOVAR(I2)
    A3=ZERO
    DO 20 IKHI=1,2
    A3=A3+TOF(IKHI,I1)*TOCOVA(IKHI,I2)
20 CONTINUE
    A4=TOCOVD(IALF,I3,ILAM)
    SOM=A1*(A2+A3+A4)
    PRINT 1100,IALF,IBET,ILAM,IMUU,A1,A2,A3,A4,SOM
4 CONTINUE
    SOM=ACOV(1)*COUR(2)+ACOV(2)*(COUR(4)-COUR(1))-ACOV(3)*COUR(3)
    TOM=ACOV(1)*COUR(4)+ACOV(2)*(COUR(6)-COUR(3))-ACOV(3)*COUR(5)
    CALL AVANPA(0)
    PRINT 1101,SOM,TOM
    RETURN
END
SUBROUTINE MIXTE(ACONTR,A,B,C,CAB)

IMPLICIT REAL*8(A-H),REAL*8(O-Z)

DIMENSION ACONTR(3),A(3),B(2,2),C(3),CAB(2,2)

DIMENSION ACT(2),T(2,2)

ACT(1,1)=ACONTR(1)

ACT(1,2)=ACONTR(2)

ACT(2,1)=ACONTR(2)

ACT(2,2)=ACONTR(3)

BCT(1,1)=A(1)

BCT(1,2)=A(2)

BCT(2,1)=A(2)

BCT(2,2)=A(3)

C CALCUL DE B...MIXTE
CALL MATA(ACT,BCT,B)

C CALCUL DE B...CONTRA
CALL MATA(ACT,B,BCT)

C CALCUL DE LA FORME CONTRACTEE
CALL MATB(B,B,CAB)

C(1)=BCT(1,1)

C(2)=BCT(2,1)

C(3)=BCT(2,2)

RETURN

END

DOUBLE PRECISION FUNCTION PROSCA(N,A,B)
DOUBLE PRECISION ZERO,A,B
DIMENSION A(N),B(N)
DATA ZERO/0.DO/

PROSCA=ZERO

DO 1 I=1,N

1 PROSCA=PROSCA+A(I)*B(I)

RETURN

END

SUBROUTINE PROVE1(X,Y,Z)
IMPLICIT REAL*8(A-H),REAL*8(O-Z)
DIMENSION X(3),Y(3),Z(3),IA(4)
DATA IA/2,3,1,2/
DO 1 I=1,3
I1=IA(I)
I2=IA(I+1)
Z(I)=X(I1)*Y(I2)-X(I2)*Y(I1)

1 CONTINUE

RETURN

ENTRY PROVE2(X,Y,Z)

DO 2 I=1,3

I1=IA(I)

I2=IA(I+1)

Z(I)=Z(I)+X(I1)*Y(I2)-X(I2)*Y(I1)

2 CONTINUE

RETURN

END
SUBROUTINE SORTIE(PHI,ACOVAR,ACONTR,ADS,BCOVAR,BCONTR,BMIXT,TOF,CA
18,COUR,AL,AM)

IMPLICIT REAL*8(A-H),REAL*8(O-Z)

COMMON/BURN/

ALPHA,BETA

DIMENSION PHI(3,2),ACOVAR(3),ACONTR(3),BCOVAR(3),BCONTR(3),IA(3,2)

DIMENSION TOF(2,3),BMIXT(4),COUR(2,3),CAB(4),AM(4,12),AL(4,12)

1100 FORMAT(1H,3PHI,2X,1HO,9X,2H 1,8X2H 2,8X,2H11,8X,2H12,8X,2H22,8X,
13H1*2,7X,4H11*1,6X,4H12*2,6X,3H111,7X,3H112,7X,3H122,7X,3H222/1H ,
23(1H*)/1H-)

1101 FORMAT(1X,13(1X,G9.2))

1102 FORMAT(10X,3HADS,1X,D9.2,23X,0HMETRIQUE,27X,5HXSI 1,3X,D15.8,3X,5H
1XSI 2,5X,D15.8/10X,3(1H*),33X,8(1H*),27X,5(1H*),21X,5(1H*))

1103 FORMAT(12(1X,D10.3))

1104 FORMAT(30X,D15.8,10X,D15.8)

1105 FORMAT(40X,37HCOURBURE COVARIANTE ET CONTRAVARIANTE/40X,37(1H*))

1106 FORMAT(50X,33HCOMPOSANTES MIXTES DE LA COURBURE/50X,33(1H*))

1107 FORMAT(10X,8(1X,D14.7))

1108 FORMAT(10X,6(2X,D15.8))

1109 FORMAT(54X,24HSYMBOLES DE CHRISTOFFELN/54X,24(1H*))

1110 FORMAT(54X,23HDERIVEES DE LA COURBURE/54X,23(1H*))

1111 FORMAT(54X,15HEDITION DES AL/54X,15(1H*)/1X)

1112 FORMAT(54X,15HEDITION DES AM/54X,15(1H*)/1X)

1113 FORMAT(50X,32HCONTROLE DE LA RELATION DE GAUSS/50X,32(1H*)/1X)

CALL AVANPA(0)

PRINT 1100

DO 50 I=1,3

PRINT 1101,{PHI(I,J),J=1,13}

50 CONTINUE

CALL AVANPA(0)

PRINT 1102,ADS,ALPHA,BETA

CALL AVANPA(0)

PRINT 1104,{ACOVAR(I),ACONTR(I),I=1,3}

CALL AVANPA(0)

PRINT 1105

CALL AVANPA(0)

CALL AVANPA(0)

PRINT 1104,{BCOVAR(I),BCONTR(I),I=1,3}

CALL AVANPA(0)

PRINT 1106

CALL AVANPA(0)

PRINT 1107,{BMIXT(I),I=1,4),(CAB(I),I=1,4)

CALL AVANPA(0)

PRINT 1109

CALL AVANPA(0)

PRINT 1108,{(TOF(I,J),J=1,3),I=1,2}

CALL AVANPA(0)

PRINT 1110

CALL AVANPA(-1)

PRINT 1108,{(COUR(I,J),J=1,3),I=1,2}

CALL AVANPA(1)

PRINT 1111

CALL NOLMAT(4,12,AL)

DO 51 I=1,4

51 PRINT 1103,{AL(I,J),J=1,12}

CALL AVANPA(0)

PRINT 1112
CALL NOLHAT(4,12,AM)

DO 52 I=1,4

52 PRINT 1103,{AM(I,J),J=1,12}

CALL AVANPA(0)

PRINT 1113

CALL METGOS(PHI,TOF,BCOVAR,ACOVAR,COUR)

RETURN

END

SUBROUTINE STOFFE (PHI,ACONTR,TOFFEL)

IMPLICIT REAL*8(A-H),REAL*8(O-Z)

DIMENSION PHI(3,7),AUX(3),ACONTR(3),TOFFEL(2,3)

DO 1 I=1,3

AUX(I)=ACONTR(1)*PHI(I,2)+ACONTR(2)*PHI(I,3)

1 CONTINUE

DO 2 I=1,3

2 TOFFEL(1,I)=PROSCA(3,AUX,PHI(1,I+3))

DO 3 I=1,3

AUX(I)=ACONTR(2)*PHI(I,2)+ACONTR(3)*PHI(I,3)

3 CONTINUE

DO 4 I=1,3

4 TOFFEL(2,I)=PROSCA(3,AUX,PHI(1,I+3))

RETURN

END