Dw : 03324 RR248

IRIA
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

<figure><img src="image_1.png" /></figure>

SOME EXAMPLES
OF IMPLEMENTATION
AND OF APPLICATION OF
THE FINITE ELEMENT METHOD

Bertrand MERCIER
Olivier PIRONNEAU

Rapport de Recherche N° 248

Août 1977

A4 9p.
SOME EXAMPLES
OF IMPLEMENTATION AND OF APPLICATION
OF THE FINITE ELEMENT METHOD

Bertrand MERCIER and Olivier PIRONNEAU

IRIA / LABORIA

Résumé :

On décrit un certain nombre de méthodes pour implémenter sur ordinateur la méthode des Éléments Finis.

A titre d'illustration, on donne des exercices et des programmes FORTRAN simples, prêts à être
utilisés par le lecteur et travaillant en mémoire centrale.

Le plan est le suivant :

1) Rappels sur l'optimisation avec ou sans contraintes.

2) 'implémentation des éléments finis les plus simples (avec gradient conjugué).

3) Cas général : triangulations, matrices de rigidité, matrices élémentaires.

4) Applications : problèmes non linéaires.

5) Applications : inéquations variationnelles.

6) Applications : contrôle optimal et optimum design.

Abstract:

We describe some methods of implementation of the finite elements method. We give some exercises and some simple FORTRAN programs ready to be used by the reader.

An outline of these notes is the following :

1) Some recalls on unconstrained and constrained optimization.

2) Implementation of the simplest finite elements (with conjugate gradient).

3) General case: triangulations, stiffness matrices, sub-matrices.

4) Applications : non linear problems.

5) Applications : varional inequalities.

6) Applications : optimal control and optimum design.
<figure><img src="image_2.png" /></figure>

## Foreword

These notes are complementary to those of LIONS. Their purpose is to help the research engineer implement some finite element programs on a medium size computer.

We have also taken this opportunity to show some applications of the finite element method.

We begin in chapter 1 with a brief presentation of the most important methods in optimization, and for the solution of simultaneous linear equations.

A detailed explanation of the programming of the simplest finite element method can be found in chapter 2.

In chapter 3 we describe the various steps for the construction of a program in more general circumstances. (Chapters 2 and 3 can be read independently)

In chapter 4 some applications are given for non-stationary or non linear partial differential equations.

In chapter 5 we show how to solve a variational inequality, and the last chapter deals with the implementation of optimal control and optimum design.

The symbol □ means a logical end of a proof, exercise, or algorithm.

To relieve the reader from unnecessary burden of programming we give a listing of some Fortran subroutines useful for finite elements. These are
CERCLE triangulation of a Ring

SQUARE Square

DPF one line minimization algorithm

FEMIA computes $\int_{\Omega} (\nabla u \cdot \nabla v + a_o uv)dx$ when u, v are piece wise linear

FEMIBF computes $\int_{\Omega} fv dx$

FEMIBG computes $\int_{\Gamma} gv d\Gamma$

FEMICG solves the mixed Neuman-Dirichlet problem for Poisson's equation discretized by the FEM of degree 1. The variational form is minimized by the conjugate gradient method

FEMICF same as above but the linear system is solved by Choleski factorization

FEM2RM solves the Dirichlet problem by the FEM of degree 2. The linear system is solved by the relaxation method

TRI Refines a triangulation

A Fortran program for the solution of an optimal control problem is given in chapter 6.

We thank Dr. Bourgat for providing FEM2RM.
# PLAN

## I. Preliminaries

1. Optimization

- The steepest descent method (SDM)

- A one line search algorithm

- The conjugate gradient method (CGM)

2. Constrained optimisation

- The penalty method

- The augmented Lagrangian method

3. Solution of linear systems

- The Choleski factorization

- Gauss-Silvestrini's relaxation method

## II. Implementation of the finite element method: 1 Lagrangian element P1

1. Statement of the problem

2. Solution by the SDM

- Algorithm

- Implementation

3. Solution by Choleski factorization

## III. Implementation of the finite element method: 2 the general method

## IV. Examples of application

1. A linear parabolic problem

2. The wave equation
3. A non linear elliptic problem: compressible potential flow

4. A linear elliptic problem: Stokes flow

5. The incompressible Navier Stokes equation

6. The biharmonic problem

# V. Resolution of a variational inequality

# VI. Optimal Control and Optimum design

1. Resolution of an optimal control problem

2. Some examples of applications of optimal control

- control of the temperature of an oven

- integration of a non linear PDE

3. Implementation of optimum design
# CHAPTER I

## PRELIMINARIES

### 1. Optimization

#### 1.1. General Remarks

In this section we shall consider the following problem:

Find $u \in H$ such that

$$ (1) \qquad J(u) = \inf_{v \in H} J(v) $$

where $H$ is a Hilbert space, with scalar product $(\cdot, \cdot)$, and norm $\|\cdot\|$; and $J: H \to R$ a real function defined on $H$ and taking its values in the real line $\mathbb{R}$.

In many practical applications $H$ is finite dimensional, and we have $H = \mathbb{R}^N$ for some positive $N$.

A very useful information about the function $J$ is the knowledge of its gradient. The function $J$ is differentiable in $u \in H$ if there exists an element $J'(u) \in H$ such that

$$ (2) \qquad (J'(u), v) = \lim_{\lambda \to 0} \frac{1}{\lambda} [J(u + \lambda v) - J(u)] \quad \forall v \in H $$

Note that this defines a mapping $J': H \to H$. If this mapping is continuous, we shall say that the function $J$ is continuously differentiable.
**Example 1:** We recall that when $H = \mathbb{R}^N$ the element $J'(u)$ can be represented by a line-vector:

$$J'(u) = \left( \frac{\partial J}{\partial x_1}, \frac{\partial J}{\partial x_2}, \dots, \frac{\partial J}{\partial x_N} \right)$$

Where the scalar quantities $\frac{\partial J}{\partial x_i}$, which depend on $u$, for $i = 1, \dots, N$, are called **partial derivatives** of $J$. $\square$

**Example 2:** Let $H = H_0^1(\Omega)$ and

$$J(v) = \frac{1}{2} \int_{\Omega} |\text{grad } v|^2 dx - \int_{\Omega} f v dx$$

We check that $(J'(v),w) = \int_{\Omega} \text{gradw} \text{gradw} dx - \int_{\Omega} f v dx$.

Note that the dual of $H_0^1(\Omega)$ is a complicated Sobolev space and that it is sufficient for our purpose to know the value of $(J'(v),w)$ for any $w \in H$. $\square$

One of the most important remark in optimization theory is that if $u$ is a solution of problem (1) then one has $J'(u) = 0$ or equivalently

$$(3) \qquad (J'(u),v) = 0 \qquad \forall v \in V.$$

Note that the variational equality (3) gives again the variational formulation of the Dirichlet problem in the case of example 2.

Conversely, if $u$ is a solution of (3) and if the function $J$ is convex, then $u$ is a solution of the optimization problem (1).
We recall that the function J is convex if (and only if)
for all u, v ∈ V and 0 < λ < 1

$$
(+) \quad J(\lambda u + (1 - \lambda)v) \le \lambda J(u) + (1 - \lambda)J(v).
$$

If the function J is **strictly convex** (that is, if (h) holds
with a strict inequality) then the problem (1) has **at most one**
solution.

1.2. Steepest descent method

As an example of an optimization method for solving problem
(1) without using the gradient of J, we can mention the relaxation
method (LIONS [3, Ch. 3, §7.2]). We shall consider in the
following 2 methods using the gradient of J: the steepest descent
method and the conjugate gradient method.

Both of them are iterative and consist of solving at each
step a simpler problem which is the minimization of a function in
one variable, which is the restriction of the function J to a
special direction.

In the steepest descent method, the direction chosen at
each step (for the minimization in one variable) is the direction
of the gradient. This gives the following algorithm:

<caption>**Algorithm 1: Steepest descent method (S.D.M.)**</caption>

$u^0 \in H$ being given; by induction

$u^n \in H$ being given.

a. Compute $J'(u^n)$.
b. Compute $\lambda^n$ solution of

$$ \text{Inf}_{\lambda \in \mathbb{R}} \phi(\lambda) \quad \text{where} \quad \phi(\lambda) = J(u^n - \lambda J'(u^n)). $$

c. Set $u^{n+1} = u^n - \lambda^n J'(u^n)$. $\square$

We shall make the assumption that the function J is _elliptic_,
that is, for some positive $\alpha$,

$$ (E) \qquad (J'(u) - J'(v), u - v) \ge \alpha \|u - v\|^2 \quad \forall u, v \in V. $$

As a matter of fact by using the identity

$$ \psi(1) - \psi(0) = \int_{0}^{1} \psi'(t) dt \quad \text{where } \psi(t) = J(u + t(v - u)), \text{ one can} \\ \text{prove that (E) implies} $$

(*) $J(v) \ge J(u) + (J'(u), v - u) + \frac{\alpha}{2} \|u - v\|^2$, $\forall u, v \in H$, which proves that the problem (1) has a unique solution that we shall denote by $u$ (see LIONS [3, Ch. 3, Remark 1.2]).

We have the following convergence result:

**Theorem 1:** Under the assumption (E), if J' is lipschitz continuous (i.e., for some c, $||J'(v) - J'(w)|| \le c||v - w||$, $\forall v, w$), then the steepest descent algorithm converges:

$$ u^n \to u \text{ in } H(\text{strongly}) \text{ as } n \to +\infty. $$

**Proof:** The scalar $\lambda^n$ is characterized by the equation

$$ \phi'(\lambda^n) = (J'(u^n - \lambda^n J'(u^n)), J'(u^n)) = (J'(u^{n+1}), J'(u^n)) = 0. $$
As a matter of fact, this proves that _two successive descent directions are orthogonal_ in the steepest descent algorithm.

Hence

$$ (J'(u^{n+1}), u^{n+1} - u^n) = 0 $$

so that, applying (*), we get

$$ J(u^n) - J(u^{n+1}) \ge \frac{\alpha}{2} \| u^n - u^{n+1} \|^2 . $$

As the sequence $J(u^n)$ is decreasing and bounded from below, it is converging, hence $\| u^n - u^{n+1} \| \to 0$, as $n \to \infty$.

Furthermore, applying the orthogonality of $J'(u^n)$ and $J'(u^{n+1})$:

$$ \| J'(u^n) \|^2 = (J'(u^n), J'(u^n) - J'(u^{n+1})) $$

$$ < C \| J'(u^n) \| \| u^n - u^{n+1} \| $$

where we applied the fact that $J'$ is Lipschitzian. Applying now $J'(u) = 0$ and (E) once again, we get, as $J'(u^n) \to 0$:

$$ \| u^n - u \| \to 0. \quad \square $$

Note that weaker assumptions are sufficient if $H = \mathbb{R}^N$ (see CEA [9]).

### Exercise 1:

Use the steepest descent method for solving problem 1 with $H = \mathbb{R}^2$ and $J(u) = J(u_1, u_2) = a_1 u_1^2 + a_2 u_2^2$

(Step b of the steepest descent algorithm can be solved explicitly). Choose first

$$ a_1 = 1, a_2 = 2 $$

$$ \text{then } a_1 = 1, a_2 = 100. \quad \square $$
**Exercise 2:** Compute the gradient of the quadratic function
$J: \mathbb{R}^N \to \mathbb{R}$ defined by

$$J(u) = \frac{1}{2} u^T A u - b^T u$$

where A is a N x N symmetric matrix, and b a column vector.

Let $w \in \mathbb{R}^N$ be a given descent direction compute explicitly
$\lambda$ minimizing $\phi(\mu) = J(u - \mu w)$ for $\mu \in \mathbb{R}$. $\square$

### 1.3. Methods for solving a minimization problem in one variable

The problem of minimizing a function $\phi$ in one variable,
which appears at step b of the steepest descent algorithm, can
be explicitly solved if $\phi$ is a polynomial of degree less than 4.

Otherwise, we can use any method for solving the equivalent
non linear equation (in one variable).

$$\phi'(\lambda) = 0.$$

Note that the same minimization problem in one variable
occurs in the relaxation method described in LIONS [3, Ch. 3, §7.2]

If it is possible to compute $\phi''$ (that is, in a practical way,
if we can compute the Hessian of the function $J$), then one can use
the Newton algorithm to solve $\phi'(\lambda) = 0$.

Otherwise a very good method is to search by trials and
errors $\lambda$ and $\mu$ such that

$$\phi\left(\frac{\lambda + \mu}{2}\right) \le \min(\phi(\lambda), \phi(\mu))$$
and to choose $\lambda^n$ achieving the minimum of the parabola interpolating
the 3 points $(\lambda, \phi(\lambda))$, $(\mu, \phi(\mu))$ and $(\frac{\lambda+\mu}{2}, \phi(\frac{\lambda+\mu}{2}))$. (See Fortran
Subprogram DPF Further.)

To search for $\lambda$ and $\mu$, one can use the following
**Algorithm** (Dichotomy)

$\lambda^0, \mu^0 \in \mathbb{R}$ being given; by induction

$\lambda^h, \mu^h \in \mathbb{R}$ being given.

$$
a. \quad \underline{\text{Compute}} \quad \alpha = \phi(\lambda^h), \beta = \phi(H^h), \gamma = \phi\left(\frac{\lambda^h + \mu^h}{2}\right)
$$

$$
b. \quad \underline{\text{If}} \ \alpha > \gamma > \beta \quad \underline{\text{set}} \ \lambda^{h+1} = \frac{\mu^h + \lambda^h}{2}, \ \mu^{h+1} = \frac{\mu^h + (\mu^h - \lambda^h)}{2}
$$

$$
\underline{\text{If}} \quad \alpha < \gamma < \beta \quad \underline{\text{set}} \quad \lambda^{h+1} = \frac{\lambda^h - (\mu^h - \lambda^h)}{2}, \quad \mu^{h+1} = \frac{(\lambda^h + \mu^h)}{2}
$$

$$
\underline{\text{If}} \quad \gamma > \alpha > \beta \quad \underline{\text{set}} \quad \lambda^{h+1} = \frac{\lambda^h}{2}, \mu^{h+1} = \frac{(\lambda^h + \mu^h)}{2} \quad \underline{\text{and go back to step }} \alpha

\underline{\text{Otherwise}} \quad \underline{\text{stop.}} \quad \square
$$

1.4. **<u>Conjugate gradient method:</u>**

In the case of a quadratic functional (see exercise 2 above),
the convergence speed can be greatly increased by choosing another
descent direction than the gradient.

As a matter of fact, the sequence of descent directions
will be construct as so called conjugate directions (conjugate
for the matrix A in the case of exercise 2), in such a way that
the method converges in at most N steps for a quadratic function
in R<sup>N</sup>.

The algorithm is the following
**Algorithm 2:** conjugate gradient method (C.G.M.)

$u^0$ H being given, set $w^0 = r^0 = J'(u^0)$. By induction

$u^n, w^n$ H being given,

a. compute $r^{n+1} = -J'(u^n)$

b. compute $w^{n+1} = r^{n+1} + \theta^n w^n$ where $\theta^n = \frac{(r^{n+1} - r^n, r^{n+1})}{(r^n, r^n)}$

c. compute $\lambda^n$ solution of

$$ \underset{\lambda \in \mathbb{R}}{\operatorname{Inf}} J(u^n - \lambda w^{n+1}) $$

d. Set $u^{n+1} = u^n - \lambda^n w^{n+1}$. □

As a matter of fact, the conjugate gradient method may converge quite faster than the steepest descent method even when the function $J$ is not quadratic. However in such a case the method may fail to converge in less than $N$ iterations, and it is necessary to set $w^{N+1} = r^{N+1}$, to start again in the gradient direction.

Note that the second method necessitates very few more computations than the first one.

**Exercise 3:** Same exercise than exercise 1 but using the conjugate gradient algorithm. □

## 2. Constrained optimization

We shall consider problems of the following type:

$$ (\mathcal{P}) \qquad \underset{v \in K}{\operatorname{Inf}} J(v) $$
Where $J:R^N \to R$ is the so-called cost function, which we
shall suppose convex and continuously differentiable, and $KCR^N$
is a closed convex subset of $R^N$.

We shall denote by $\| \cdot \|$ the euclidean norm of $R^N$, and
by (,) the associate scalar product. At last:

$\pi:R^N \to K$

is the projection operator onto $K$: defined by

$$
(5) \quad (u - \pi u, \pi u - v) \ge 0 \quad \forall v \in K
$$

2.1. The gradient algorithm with projection

This is the following algorithm:

$$
\text{(6)} \quad \left\{
\begin{array}{@{}l@{}}
\text{o} \in K \text{ and } p > 0 \text{ being given} \\
\text{n} \in K \text{ being given by induction, compute} \\
\text{n} + 1 = \pi(\text{n} - p\text{J}'(\text{n}))
\end{array}
\right.
$$

Assume that the gradient $J'$ is lipschitz continuous, that is

$$
(7) \quad \|J'(u) - J'(v)\| \le C \|u - v\|, \quad \forall u, v \in \mathbb{R}^N
$$

Then, we have the following result due to GOLDSTEIN [2]:

Theorem 2: Assume that $\bar{c}$ has at least one solution, and that
(7) holds for some constant $c > 0$, then the algorithm (6)
converges for $0 < p < \frac{1}{2c}$ in the following sense: the sequence
{u^n} remains bounded, and there exists a subsequence converging to a solution of P.

We shall prove the theorem under the stronger assumption: there exists $\alpha > 0$ such that $\forall u, v \in V$

$$ (8) \quad (J'(u) - J'(v), u - v) \ge \alpha \| J'(u) - J'(v) \|^2. $$

(Note that (8) implies (7) with $C = \frac{1}{\alpha}$)

**Proof:** The solution $u$ of P can be characterized by the variational inequality (cf. LIONS, Ch. 3, §1):

$$ (J'(u), v - u) \ge 0 \quad \forall v \in K $$

which implies as a special case ($\forall p \in \mathbb{R}$) (apply (5))

$$ u = \pi(u - pJ'(u)) $$

Then

$$ u^{n+1} - u = \pi(u^n - pJ'(u^n)) - \pi(u - pJ'(u)). $$

Using the contraction property $|| \pi a - \pi b || < || a - b ||$, we have:

$$ || u^{n+1} - u ||^2 < || u^n - u - p(J'(u^n) - J'(u)) ||^2 $$

then

$$ || u^{n+1} - u ||^2 < || u^n - u - p(J'(u^n) - J'(u)) ||^2 $$

then
$$ \begin{gathered} || u^{n+1} - u ||^2 < || u^n - u ||^2 - 2p(J'(u^n) - J'(u)), \\ u^n - u) + p^2 || J'(u^n) - J'(u) ||^2. \end{gathered} $$

Using the assumption (8), we get

$$ || u^{n+1} - u ||^2 < || u^n - u ||^2 + p(p - 2\alpha) || J'(u^n) - J'(u) ||^2 $$

This proves that the sequence $\{|| u^n - u ||\}$ is decreasing for

$0 < p < 2\alpha$. The sequence $\{u^n\}$ remaining bounded we can extract a subsequence converging to $w \in K$. As we have $w = \pi(w - pJ'(w))$ and then $w$ is a solution of $F$. □

**Remark 1:** Note that the proof is still true if we change the scalar product (and the form of $R^N$). Let

$$ (u, v) = u^T A v \quad (\text{instead of } u^T v) $$

where $A$ is a symmetric and positive definite matrix.

The main effect of this change is to replace the gradient $J'(v)$ by $A^{-1}J'(v)$. The new method is still a gradient method, but with an **auxiliary operator**.

This change can improve very much the convergence in some cases where the Hessian of $J$ is "not too far" from $A$. □

An important case where (8) holds is the case where the **conjugate** $F$ of $J$, defined by

$$ (10) \quad F(v) = \sup_{w \in R} [(w,v) - J(w)], $$
is coercive. An example of such a situation is the case where

$$F(v) = \frac{\alpha}{2} \|v\|^2 + j(v)$$

where $j$ is convex and lower semi continuous (1.s.c.)

**Remark 2:** From the proof of the theorem we can see that

$$J'(u^n) \to J'(u)$$

when $0 < p < 2\alpha$. In fact from (8) again, it is possible to prove that $J'(u)$, where $u$ is solution of $J'$, is unique, even though $u$ may not be unique. □

**Exercise 4:** Consider the case where

$$(J'(u) - J'(v), u - v) \ge \alpha \|u - v\|^2 \quad \forall u, v \in V$$

and adapt the proof of the theorem. □

## 2.2. The penalty method

The projection on the convex set $K$ can be a difficult problem if $K$ is complicated.

We shall consider in this section the case where

$$K = \{v \in \mathbb{R}^N : G_j(v) < 0 \quad 1 \le j \le m\}$$

and the functions $G_j : \mathbb{R}^N \to \mathbb{R}$ are convex continuous, so that

$K$ is closed.
In such a case it is difficult to compute the projection
on the convex K.

We shall consider the function $P: \mathbb{R}^N \to \mathbb{R}$ defined by

$$
P(v) = \sum_{j=1}^{m} (G_j^+(v))^2 \equiv \sum_{j=1}^{m} (\max(0, G_j(v)))^2
$$

We notice that P is convex and positive. Moreover

$$
P(v) = 0 \Leftrightarrow v \in K.
$$

The quantity P can then be viewed as a "penalty" term.

We consider now the "penalized" problem:

$$
(P_{\epsilon}) \qquad \inf_{v \in \mathbb{R}^N} J_{\epsilon}(v)
$$

where

$$
J_{\epsilon}(v) \equiv J(v) + \frac{1}{\epsilon} P(r) .
$$

**Theorem 3:** *Assume that J is strictly convex, and that*

$$
J(v) \to +\infty \quad \text{as} \quad v \to +\infty
$$

then $P_{\epsilon}$ (resp. $P$) has unique solution $u_{\epsilon}$ (resp. $u$).

Moreover $u_{\epsilon} \to u$ as $\epsilon \to 0$.

Proof:

From LIONS [3, Ch. 3, Remark 1.2] we know that $P_{\epsilon}$ and
$\epsilon$ have a unique solution. Furthermore, we have
$$
(11) \quad J(u_{\epsilon}) \in J(u_{\epsilon}) \in \inf_{v \in K} J(v) = \inf_{v \in K} J(v) = J(u)
$$

hence $u_{\epsilon}$ remains bounded in $R^{N}$, $\forall \epsilon > 0$. We can then construct
a decreasing sequence $\epsilon_n \to 0$ as $n \to \infty$ such that

$$
u_{\epsilon} \to w \quad \text{in } R^{N} \text{ as } n \to +\infty.
$$

From (11) we deduce that

$$
(12) \quad J(w) \le J(u).
$$

and tha

$$
\frac{1}{\epsilon_n} P(u_{\epsilon_n}) \le J(u) - J(u_{\epsilon_n}).
$$

Then $P(u_{\epsilon_n}) \to 0$ as $n \to +\infty$ which proves that

$$
P(w) = \lim_{n \to +\infty} P(u_{\epsilon_n}) = 0.
$$

Hence $w \in K$, which, together with (12) proves that $w = u$.

At last the uniqueness proves that the whole sequence $u_\epsilon$
converges to $u$. $\square$

Note that for $\epsilon > 0$, $u_\epsilon$ may not belong to K: the method
is then called an **exterior penalty method**.

For other types of penalty method, we refer the reader
to GLOWINSKI-LIONS-TREMOLIERES [4].
## 2.3. Duality methods

### a) a linear situation

We shall consider first a linear problem in order to
define some vocabulary: consider the problem

$$ (13) \quad A_p + B_u^T = 0 $$

$$ (14) \quad B_p = b $$

where $A$ is a $nxn$ matrix, symmetric and positive definite

$B$ a $mxn$ matrix of maximal rank ($m < n$).

$b$ a column vector with $m$ components.

Let us try to eliminate $p$: multiplying (13) by $A^{-1}$, we get:

$$ p = -A^{-1}B_u^T $$

and then we get that $u$ is solution of the linear system

$$ (15) \quad B A^{-1} B_u^T = -b $$

which has a unique solution since the matrix $BA^{-1}B^T$ is positive
definite ($v \neq 0$ implies $v^TBA^{-1}B^T = y^T A^{-1}y$ which is positive
since $y = B^T v \neq 0$).

It is a little bit more difficult to eliminate $u$: let

$$ E = \{ \eta \in R^n : B_i = b \}. $$
From (13) we get for all $q \in E$:

$$ (Ap, q - p) \equiv (q - p)^T Ap = -(B(q - p))^T u = 0 $$

Therefore $p$ minimizes $\frac{1}{2}(Aq, q)$ on $E$, which we can write under the form of an optimization problem: $p$ is solution of

$$ (17) \qquad \text{Inf}_{q \in E} \frac{1}{2}(Aq,q). $$

If we have then an initial problem (or primal problem) with constraints of the form (14), it may happen that its dual problem (here (15)) be easier to solve.

Here, the dual problem (15) is obviously easier to solve since it has no constraints, and is linear.

We remark that the primal and the dual problem have no apparent connection. However both are a special form of the initial system (which is in some sense the real problem). The equation (13) relating the primal variable $p$, to the dual one $u$, is called extremality relation, and has a particular importance.

b) Non quadratic function with linear constraints

We can now consider a generalization of the problem (17) to a non linear case:

$$ (18) \qquad \text{Inf}_{q \in E} J(q) $$

Where $J: \mathbb{R}^N \to \mathbb{R}$ is a convex (non quadratic) function, and $E$ is still defined as in (16) (Then the constraints are still linear)
We shall try to find the dual problem of (18).
We notice first that

$$
\sup_{v \in \mathbb{R}^m} [-v^T B_q + v^T b] = \begin{cases} 0 & \text{if } q \in E \\ +\infty & \text{otherwise;} \end{cases}
$$

then (18) is equivalent to

$$
\inf_{v \in \mathbb{R}^n} \sup_{q \in \mathbb{R}^m} L(q,v)
$$

where $L: \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}$ is the following function (called Lagrangian):

$$
L(q,v) = J(q) - v^T B_q + v^T b .
$$

A general way of defining a dual problem to a given problem
is to achieve the permutation of the infinum and the suprenum.

We get:

$$
(19) \quad \sup_{v \in \mathbb{R}^m} \inf_{q \in \mathbb{R}^n} L(q,v)
$$

Let $F$ be the conjugate of $J$ defined as in (10), we see
that the dual problem (19) is equivalent (by changing signs)
to find $u$ solution $u^*$

$$
(20) \quad \inf_{v \in \mathbb{R}^m} [F(B_v^T) - b^T v] .
$$

We can check that (20) is equivalent to (15) when

$$
J(q) = \frac{1}{2} q^T A q
$$
In fact, $z \in \mathbb{R}^n$, we have

$$F(z) = \sup_{q \in \mathbb{R}^n} [(z,q) - J(q)] = \sup_{q \in \mathbb{R}^n} [z^T q - \frac{1}{2} q^T A q].$$

The last supremum is reached for $q = A^{-1}z$ and its value is

$$F(z) = \frac{1}{2} z^T A^{-1} z,$$

which sets that (15) is a special case of (20).

The numerical experiments have shown that it is more interesting to try to solve the primal problem (17) if J is differentiable, and the dual problem (20) if the conjugate F is differentiable.

Of course the primal problem has the constraint $q \in E$, but it is easy to compute the projection on E and therefore to use the *gradient algorithm with projection* for solving (17).

In fact, let $P = B^T (BB^T)^{-1} B$ denote the projection on the range of $B^T$, the projection $\pi$ from $\mathbb{R}^n$ onto E is given by

$$\pi q = (I - P)q + B^T (BB^T)^{-1} b.$$

**Remarks:** 1. If you can prove that one of the problems (17) or (20) has a unique solution, the other is differentiable! (See EKELAND-TEMAM) [5]

2. The extremality relation between u and p is now:

$J'(p) = B^T u$ if J is differentiable

or $F'(B^T u) = p$ if F is differentiable.
Note that each of these relations is a generalization
of the relation (13), which was the extremality relation in
the linear case. □

**Examples:** The previous concepts can be extended to infinite dimensional situations: replace $R^n$ and $R^m$ by 2 Hilbert spaces $Y$ and $V$, $A \in L(Y,Y')$; $B \in L(Y,V')$ (see EKELAND-TEMAM). The validity of the exchange of the infimum and the supremum is more difficult to justify. We consider the case where $Y = (L^2(\Omega))^2 = Y'$; $V = H_0^1(\Omega)$; $A$ is the identity and $B$ the operator gradient, and $B$ the divergence.

The equivalent of the system (13)(14) is

$$
(21) \quad \mathbf{p} - \operatorname{grad} \mathbf{v} = \mathbf{0}
$$

$$
(22) \quad \operatorname{div} \mathbf{p} + \mathbf{b} = \mathbf{0}
$$

where $b$ is given in $V'$.

In such a case the dual problem (15) is the dirichlet
problem: find $u \in H_0^1(\Omega)$ such that

$$
-\Delta u = b .
$$

The primal problem is

$$
\text{Inf}_{q \in E} \quad \frac{1}{2} \|q\|^2
$$

where

$$
(23) \quad E = \{ \mathbf{p} \in (L^2(\Omega))^2 : \operatorname{div} \mathbf{p} + \mathbf{b} = \mathbf{0} \}
$$
and $||p||^2 = \int_{\Omega} |p(x)|^2 dx$ ($|\cdot|$ denotes the extension norm of $R^2$)

We shall consider 2 examples of non linear cases:

**Example 1: visco-plastic problem (The flow of a Bingham fluid)**

We define (for $g > 0$)

$$C = \{q \in L^2(\Omega)^2: |q(x)| < g \text{ a.e.}\}$$

and $\pi:(L^2(\Omega))^2 \rightarrow C$ the projection onto C defined almost
everywhere as follows

$$\pi q(x) = \begin{cases} q(x) & \text{if } |q(x)| < g. \\ g \frac{q(x)}{|q(x)|} & \text{otherwise.} \end{cases}$$

We consider then the problem

$$ (24) \quad \begin{array}{ll} \text{Inf} & J(q) \\ \qquad \text{q} \in E & \end{array} $$

where $J(q) = \frac{1}{2} ||q - \pi q||^2$, and $E$ is as in (23).

We can prove that the conjugate of $J$ is $F$, defined for

$z \in (L^2(\Omega))^2$ by

$$ (25) \quad F(z) = \frac{1}{2} ||z||^2 + j(z) $$

where

$$ j(z) = g \int_{\Omega} |z(x)| dx. $$
The dual problem is then

$$ \inf_{v \in H_0^1(\Omega)} \left[ \frac{1}{2} \| \text{grad } v \|^2 + j(\text{grad } v) - \langle b, v \rangle \right] $$

where <, > is duality between V' and V.

Note that this problem is exactly the one studied by LIONS [3, Ch. 3, 8.1].

In this case as we noted before it will be easier to try to solve the primal problem (24) because J is differentiable:

$$ J'(q) = q - \pi q. \quad \square $$

### Example 2: Elasto-plastic problem.

We consider the case where

$$ J(q) = \frac{1}{2} \|q\|^2 + \psi_c(q) $$

where $\psi_c$ is the indicator function of C defined by

$$ \psi_c(q) = \begin{cases} 0 & \text{if } q \in C \\ +\infty & \text{otherwise.} \end{cases} $$

We still want to solve the problem

$$ \begin{array}{ll} \inf_{q \in E} & J(q) \end{array} $$
which is in fact equivalent to the problem

$$
(26) \quad \inf_{q \in C^1(E)} \frac{1}{2} \| q \|^2.
$$

This problem is a model for an elasto-plastic problem,
because you are looking for a stress field q in static
equilibrium (q∈E), plastically admissible (q ∈ C) and
minimizing an energy functional.

In this case, the function J is not differentiable. As
(26) has a unique solution, you can expect F to be differentiable

In fact

$$
F(z) = \int_{\Omega} \phi(z(x)) \, dx
$$

where $\phi: \mathbb{R}^2 \rightarrow \mathbb{R}$ is defined as follows:

$$
\phi(z) = \begin{cases} \frac{1}{2} |z|^2 & \text{if } |z| < \delta \\ g(|z| - \frac{\delta}{2}) & \text{otherwise.} \end{cases}
$$

(You can check that $\phi$ and then $F$ is differentiable).

In such a case, it will be easier to solve the dual problem
(20) than the primal problem (26). □

Of course, from a practical point of view, only the discrete
problem is solved (on the computer) and we are then reduced to
a finite dimensional case.
### c) Problems with inequality constraints

We shall consider here the application of duality methods to solve a problem with inequality constraints as in §2.2:

$$ (\circ) \quad \begin{array}{ll} \text{Inf} & J(v) \\ v \in K & \end{array} $$

$$ \text{where } K = \{ v \in \mathbb{R}^N : G_j(v) \le 0, \ 1 \le j \le m \} $$

and $G_j : \mathbb{R}^N \to \mathbb{R}$ are convex and continuous for $1 \le j \le m$.

We notice again that the quantity

$$ \sup_{\mu \in \mathbb{R}_+^m} \left[ \sum_{j=1}^m \mu_j G_j(v) \right] = \begin{cases} 0 & \text{if } v \in K \\ +\infty & \text{otherwise} \end{cases} $$

(Note that $\lambda \in \mathbb{R}_+^m \stackrel{\lambda_j > 0}{\Longleftrightarrow} \lambda_j \neq 0$. The initial problem $\mathcal{P}$ is then equivalent to)

$$ \begin{array}{lll} \text{Inf} & \text{Sup}_{\mu \in \mathbb{R}_+^m} & L(v, \mu) \\ v \in \mathbb{R}^N & & \end{array} $$

where $L$ is the following function called Lagrangian:

$$ L(v, \mu) = J(v) + \sum_{j=1}^{m} \mu_j G_j(v) . $$

To compute the solution of the problem $\mathcal{P}$, one can use the following Uzawa's algorithm (in a similar manner to LIONS [3, Ch. 3, §8.2]):

* $p > 0$ and $\lambda^0 \in \mathbb{R}_+^m$ being given

* by induction $\lambda^n \in \mathbb{R}_+^m$ being given, compute
$v^n$ **minimizing** $L(v, \lambda^n)$ **on** $R^N$

. compute $\lambda_j^{n+1} = \max(0, \lambda_j^n + pG_j(v^n))$.

The convergence of this algorithm is proved for $p$ sufficiently small under reasonable hypotheses (see ARROW-HURWICZ-UZAWA [6]). We refer the reader to ROCKAFELLAR [7] and EKELAND-TEMAM [5] for additional details, on these methods.

## §3 Solution of simultaneous linear equations

### 3.1. Introduction

Consider a system of simultaneous linear equations

$$ (27) \quad Ax = b $$

where $A$ is a $nxn$ matrix, $b$ and $x$ are column vectors; $A$ and $b$ are given, and we want to compute $x$.

When $A$ is positive definite, we know that solving (27) is equivalent to find $x$ minimizing

$$ (28) \quad J(y) = \frac{1}{2} y^T A y - b^T y $$

for $y \in R^n$.

Then, to solve (27) we can use all the iterative methods described in previous §1.

But all these methods are very much sensitive to what is called an **ill-conditioning effect**. Let $|| \cdot ||$ denote the **Euclidean norm** of $R^n$ and $||A||$ the associated matrix norm:
$$ \| A \| = \sup_{\| x \|_1 \le 1} \| Ax \|, $$

We call condition number of a matrix A, the number

$$ \mathrm{cond}(A) = \|A\| \|A^{-1}\|. $$

With such a definition it is easy to prove (and we advise the reader to do it) that, if $A \in R^N$ and

$$ A(x + \delta x) = b + \delta b $$

then one has

$$ (29) \qquad \| \delta x \| \le \mathrm{cond}(A) \| \delta b \| $$

In (29) the equality can be reached for special choices of $\| \delta b \|$. For a symmetric positive definite matrix, one has

$$ \mathrm{cond}(A) = \frac{\lambda_{\max}}{\lambda_{\min}} $$

where $\lambda_{max}$ and $\lambda_{min}$ are the largest and the smallest eigenvalue of the matrix A.

The matrices coming from finite element approximation of problems of order 2m have a condition number of order $h^{-2m}$.

This means that the condition number increases very much with the number of variables.

The number of iterations of most of the iterative methods
described previously is proportional to the condition number of the matrix A, if we apply them to the function J defined in (28).

An exception is the conjugate gradient method which converges (theoretically) in less than N iterations.

The conjugate gradient method is then a very good method which we recommend.

However, when one has to solve many times a system of equation with the same matrix A, the Gaussian elimination that we shall describe in the following is more convenient.

In fact, the result of Gaussian elimination is to write the matrix A as a product of 2 triangular matrices:

$$A = LU$$

where L is Lower triangular and U is Upper triangular.

When such a decomposition has been computed, it is very easy to solve the linear system (27) by proceeding as follows: solve first

$$(3) \quad Ly = b,$$

and then

$$Ux = y.$$

In fact the solution of a system of simultaneous linear equations, the matrix of which is triangular is very easy.
Consider the first one $L = b$ that we can write:

$$ \ell_{11} y_1 = b_1 $$

$$ \ell_{21} y_1 + \ell_{22} y_2 = b_2 $$

$$ \ell_{n1} y_1 + \ell_{n2} y_2 + \dots + \ell_{nn} y_n = b_n $$

From the first equation you deduce

$$ y_1 = \frac{b_1}{\ell_{11}} $$

$(\ell_{11}$ is different of zero, otherwise L and then A is singular).

From the second equation you deduce

$$ y_2 = \frac{1}{\ell_{22}} (b_2 - y_1) $$

which you can compute, since you know already $y_1$. By induction
at the last step, you know already $y_1$ ................ $y_{n-1}$ and you
compute $y_n$ from the last equation by

$$ y_n = \frac{1}{\ell_{n1}} (b_n - \ell_{n1} y_1 - \ell_{n2} y_2 - \dots - \ell_{nn-1} y_{n-1}) $$

## 3.2. The standard Gaussian elimination

We shall describe the standard Gaussian elimination for
factorizing A in a product of two triangular matrices:

$$ A = LU. $$
Let $A^{(1)} = A$; we shall compute by induction $A^{(k)}$ for $k = 2, \dots, m$, where

$$ (32) \qquad A^{(k)} = \left( \begin{array}{cc|c} a_{11}^{(k)} & & a_{1n}^{(k)} \\ & \ddots & \\ a_{kk}^{(k)} & & a_{kn}^{(k)} \\ \hline a_{nk}^{(k)} & & a_{nn}^{(k)} \end{array} \right) $$

We can consider that $A^{(k)}$ is the matrix of the system of simultaneous linear equations at the $k$-th step. To obtain the matrix $A^{(k+1)}$, you perform in the corresponding system of linear equations the elimination of the $k$-th unknown by adding the $k$-th equation multiplied by the coefficient

$$ - \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}} $$

to the equations $i = k+1, \dots, n$.

The $k$ first lines and the $k-1$ first columns of $A^{(k+1)}$ are chosen equal to the corresponding lines and columns of $A^{(k)}$. Furthermore

$$ (33) \qquad a_{ij}^{(k+1)} = a_{ij}^{(k)} - \frac{a_{ik}^{(k)} a_{kj}^{(k)}}{a_{kk}^{(k)}} $$

for $i = k+1, \dots, n$

$j = k, \dots, n$

provided that $a_{kk}^{(k)}$ is different of zero.
We can check that the elements $a_{ik}^{(k+1)}$ = 0 for i = k+1, ..., r so that by induction $A^{(n)}$ = U is upper triangular.

It may happen for a general non singular nxn matrix
A that $a_{kk}^{(k)}$ = 0, and if it is so, we must achieve a permutation
of the rows of A$^{(k)}$ in order to continue the Gaussian elimination.

In such a case it is not possible to write the matrix A
as a product LU, but only the matrix PA, where P is a
permutation matrix.

However this does not happen for symmetric positive
definite matrices (see WILKINSON [8]), for which all the "pivots"
$a_{kk}^{(k)}$ are different of zero.

Let now L be the lower triangular matrix

$$
L = ((\ell_{ij}))_{1 \le i j \le n}
$$

such that

$$
\ell_{ij} = 1 \quad \text{for } i = 1, \dots, n
$$

$$
\ell_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}} \quad \text{for } k = 1, \dots, n \\
\phantom{\ell_{ik} = \frac{a_{ik}^{(k)}}{}} \quad \text{and } i = k+1, \dots, n,
$$

then, you can check that A = LU, and the factorization process
is completely described.
**Remark: Cholesky Method.** As A is positive definite, we

have $a_{kk}^{(m)} > 0$ (see WILKINSON [8])

Let D be the diagonal matrix such that

$$d_{kk} = \sqrt{a_{kk}^{(m)}}$$

As A is symmetric we have

$$D^{-1}U = (DL)^T$$

Let M = DL, we have then found a lower triangular matrix M such that

$$A = MM^T. \qquad (34)$$

(Conversely if there exists a lower triangular matrix M such that (34) holds, then A is positive definite.)

The Cholesky method consists in computing directly the matrix M which is uniquely determined from the recursive formulae.

$$ (35) \qquad \begin{cases} m_{il} = a_{il} & \text{for } i = 2, \dots, n \\ m_{il} = \frac{a_{il}}{m_{li}} & \text{for } j = 2, \dots, n \end{cases} $$

for j = 2, ..., n:
$$
(36) \quad \left\{
\begin{array}{l}
m_{jj} = (a_{jj} - \sum_{k=1}^{j-1} m_{jk}^2)^{\frac{1}{2}} \\[2ex]
m_{ij} = (a_{ij} - \sum_{k=1}^{j-1} m_{ik} m_{jk}) \frac{1}{m_{jj}} \quad \text{for } i = j+1, \dots, n. \quad \square
\end{array}
\right.
$$

Remark 5: Band-Structure

For a full nxn matrix A, the Gaussian elimination requires
asymptotically $\frac{n^3}{3}$ operations (no other method can be faster in
such a case).

However, if the matrix A is sparse, considerable computing
time can be saved by avoiding unnecessary operations on zero
elements.

The process is very easy when the matrix has a "band-structure",
i.e., for some number w, called bandwidth, one has the property

$$
|i - j| > w \Rightarrow a_{ij} = 0.
$$

This kind of matrices can be represented in the following way:

<figure><img src="image_3.png" /></figure>
For such matrices the number of operations is now about
$n \frac{w^2}{2}$ for the factorization.

The number of operations is also reduced for the solution
of (30) (31). If the lower part of L is full, the solution of
$Ly = b$ requires about $\frac{n^2}{2}$ operations. For a band matrix this
number is reduced to $n$ times $w$. $\square$

**Exercise 4:** Show by using the formulae (33) that the triangular factors L, U have the same bandwidth w than A, if A has a band structure. $\square$

**Exercise 5:** Same exercise as the previous one, but for M, by using the formulae (35) (36). $\square$

**Exercise 6:** Compute the LU decomposition of the nxn tridiagonal matrix.

<figure><img src="image_4.png" /></figure>

**Exercise 7:** Show that the vector $y^\alpha \in R^n$, the components of which are given by

$$ y_i^\alpha = \sin \left( \alpha \frac{i}{n+1} \pi \right) \quad 1 \le i \le m $$

is eigenvector of the tridiagonal matrix A defined in the previous exercise, for $\alpha = 1, \dots, n$, associated to the eigenvalue
$$ \lambda_{\alpha} = 4 \sin^2 \frac{\alpha\pi}{2(n+1)} . $$

Deduce from this result that $\text{cond}(A) \sim \frac{4(n+1)^2}{\pi^2}$

for large $n$. □
REFERENCES

[1] POLAK, E.: *Computational methods in optimization*, Academic Press, 1971.

[2] GOLDSTEIN, A.A., Convex programming in Hilbert spaces, Bull. Amer. Math. Soc., Vol. 70 (1974) 709 - 710.

[3] LIONS, J.L.: These lecture notes.

[4] GLOWINSKI, R., LIONS, J.L., TREMOLIERES, R.: *Analyse Numerique Des Équations Variationnelles*, Donod, Paris, 1976.

[5] EKELAND, I., TEMAM, R.: *Analyse Convexe et Problèmes Variationnels*, Dunod, Paris, 1973.

[6] ARROW, K.J., HURWICZ, L., UZAWA, H.: *Studies in Linear and Non-linear Programming*, Stanford University Press, 1968.

[7] ROCKAFELLAR, T.: *Convex Analysis*, Princeton, 1970.

[8] WILKINSON, J.H.: *The Algebraic Eigenvalue Problem*, Clarendon Press, Oxford, 1965.

[9] CEA, J.: *Optimisation, Théorie et Algorithmes*, Dunod, Paris, 1970.
# Annex of Chapter 1

## Description of the subroutine DPF

The subroutine DPF(RO, C) computes an approximation of the minimum of $f(\rho)$.

RO and C are given on input and the subroutines looks

for an interval $\rho_0 \le \rho_1$ such that $f\left(\frac{\rho_0 + \rho_1}{2}\right) < \min(f(\rho_0), f(\rho_1))$;

then it takes the minimum of the parabolic interpolation through $\rho_0, \rho_1$, $(\rho_0 + \rho_1)/2$. The initial interval is $\rho_0 = \text{RO}/2, \rho_1 = \text{RO}$.
If $\rho_1$ becomes less than $10^{-10}$ the subroutines return.

When an evaluation of $f(\rho)$ is not costly the user may wish to have
a better precision; in that case DPF can be called several times, with
the output of the previous call as input to the next call.

DPF must be provided with a subroutine CPF (RO,F)

which computes $F = F(RO)$ for an input RO.

A listing of DPF is in annex.
ANNEX TO CHAPTER 1

Listing of DPF

RELEASE 2.0

DPF

DATE = 77060 15/33/21

PAGE 0001

SUBROUTINE DPF (RO, C)

C COMPUTES THE MINIMUM OF A FUNCTION F OF ONE VARIABLE BY DICHOTOMY

CC PLUS PARABOLIC INTERPOLATION. THE SUBROUTINE LOOKS FOR AN INTERVAL R01 R02

C SUCH THAT F ((R01 + R02)/2) LESS THAN MIN(F(R01), F(R02))

C THEN IT TAKES THE MINIMUM OF THE PARABOLIC INTERPOLATION THROUGH THE

C LAST 3 POINTS.

C INPUT RO INITIAL VALUE FOR THE DICHOTOMIC SEARCH.

C C VALUE OF THE FUNCTION AT RO=O.

C CPS (RO,F) SUBROUTINE WHICH RETURNS F=F(RO) FOR AN INPUT RO.

C OUTPUT RO APPROXIMATION OF THE OPTIMUM.

C

C CALCUL DU RO OPTIMAL PAR DICHOTOMIE

DIMENSION ROMIN(3), FMIN(3)
IS=2

<table><thead><tr><th></th><th></th><th>PORTAN IV G1</th><th>RELEASE 2.0</th><th>DPF</th><th>DATE</th></tr></thead><tbody><tr><td></td><td>XRO=1.D-10</td><td>0035</td><td>GOTO 360</td><td></td><td></td></tr><tr><td>240</td><td>CONTINUE</td><td>0036</td><td>420 CONTINUE</td><td></td><td></td></tr><tr><td></td><td>ROMIN (1)=RO*0.5</td><td>C037</td><td>425 ROMIN (1)=ROMIN (2)</td><td></td><td></td></tr><tr><td></td><td>ROMIN (2)=RO</td><td>0038</td><td>FMIN (1)=FMIN (2)</td><td></td><td></td></tr><tr><td></td><td>ROMIN (3)=RO*2.</td><td>0039</td><td>ROMIN (2)=ROMIN (3)</td><td></td><td></td></tr><tr><td></td><td>L=O</td><td>0040</td><td>FMIN (2)=FMIN (3)</td><td></td><td></td></tr><tr><td>300</td><td>L=L + 1</td><td>0041</td><td>ROMIN (3)=ROMIN (3)*2.</td><td></td><td></td></tr><tr><td></td><td>CALL CPS (ROMIN (L), FMIN (L))</td><td>0042</td><td>CALL CPS (ROMIN (3))</td><td></td><td></td></tr><tr><td></td><td>IF (IS.EQ.2) PRINT 100, ROMIN (L), FMIN (L)</td><td>0043</td><td>IF (IS.EQ.2) PRINT 100, ROMIN (3), FMIN (3)</td><td></td><td></td></tr><tr><td>100</td><td>FORMAT (1X,'POINT=', E10.3,') CRI= ',E11.4)</td><td>0044</td><td>GOTO 370</td><td></td><td></td></tr><tr><td></td><td>IF (L.NE.2) GOTO 352</td><td>0045</td><td>450 RO=ROMIN (2)</td><td></td><td></td></tr><tr><td></td><td>IF (FMIN (1).LT.FMIN(2)) GOTO 380</td><td>0046</td><td>IF (FMIN (2).EQ.FMIN (3)) GO TO 500</td><td></td><td></td></tr><tr><td>352</td><td>CONTINUE</td><td>0047</td><td>451 CONTINUE</td><td></td><td></td></tr><tr><td></td><td>IF (L.NE.1) GOTO 355</td><td>0048</td><td>C...... CALCUL DE RO INTERPOLE</td><td></td><td></td></tr><tr><td></td><td>IF (FMIN (1).LT.ECH) GOTO 355</td><td>0049</td><td>SN=0.</td><td></td><td></td></tr><tr><td></td><td>RO=RO*0.5</td><td>0050</td><td>SD=0.</td><td></td><td></td></tr><tr><td></td><td>...... TEST D'ARRET</td><td>0051</td><td>DO 454 I=1,3,1</td><td></td><td></td></tr><tr><td></td><td>IF (RO.LT.XRO) GO TO 500</td><td>0052</td><td>S=0.</td><td></td><td></td></tr><tr><td></td><td>C......</td><td>0053</td><td>PR=1.</td><td></td><td></td></tr><tr><td></td><td>GOTO 240</td><td>0054</td><td>DO 452 J=1,3</td><td></td><td></td></tr><tr><td>355</td><td>CONTINUE</td><td>0055</td><td>IF (I.EQ.J) GOTO 452</td><td></td><td></td></tr><tr><td></td><td>IF (L.LT.3) GOTO 300</td><td>0056</td><td>S=S+ROMIN (J)</td><td></td><td></td></tr><tr><td>160</td><td>IF (FMIN (1).LT.FMIN(2)) GOTO 380</td><td>0057</td><td>PR=PR*(ROMIN(I)-ROMIN(J))</td><td></td><td></td></tr><tr><td>370</td><td>IF (FMIN(2)=FMIN(3)) 450,450,420</td><td>0058</td><td>CONTINUE</td><td></td><td></td></tr><tr><td>380</td><td>ROMIN(3)=ROMIN(2)</td><td>0059</td><td>SN=SN+FMIN(I)*S/PR</td><td></td><td></td></tr><tr><td></td><td>F=IN(3)=FMIN(2)</td><td>0060</td><td>SD=SD+FMIN(I)/PR</td><td></td><td></td></tr><tr><td></td><td>ROMIN(2)= ROMIN(1)</td><td>461</td><td>RO=SN/SD/2.0</td><td></td><td></td></tr><tr><td></td><td>F=IN(2)=FMIN(1)</td><td>462</td><td>CONTINUE</td><td></td><td></td></tr><tr><td></td><td>ROMIN(1)=ROMIN(1)*0.5</td><td>468</td><td>CONTINUE</td><td></td><td></td></tr><tr><td></td><td>CALL CPS(FOMIN(1),FMIN(1))</td><td>500</td><td>CONTINUE</td><td></td><td></td></tr><tr><td></td><td>IF (IS.EQ.2)PRINT 100,ROMIN(1),FMIN(1)</td><td>501</td><td>RETURN</td><td></td><td></td></tr><tr><td></td><td></td><td>502</td><td>RETURN</td><td></td><td></td></tr><tr><td></td><td></td><td>503</td><td>RETURN</td><td></td><td></td></tr><tr><td></td><td></td><td>504</td><td>RETURN</td><td></td><td></td></tr></tbody></table>
## Chapter 2

# Implementation of the Finite Element Method I
## Lagrangian element PI

### I. Statement of the Problem

Consider the Partial Differential Equation (P.D.E.) ($a_o \ge 0$)

$$ (I) \quad \left\{ \begin{array}{l} - \Delta u + a_o u = f \quad \text{in } \Omega \\ u|_{\Gamma_-} = u_o \\ \overline{\Gamma_1 u \Gamma_2} = \partial \Omega \\ \\ \frac{\partial u}{\partial n}|_{\Gamma_2} = g \quad \Gamma_1 n \Gamma_2 = \phi \end{array} \right. $$

<figure><img src="image_5.png" /></figure>

It is shown in Lions [I] that (I) is in some sense equivalent to

$$ (2) \quad \min_{u \in U(u_o)} \left\{ \frac{1}{2} a(u, u) - b(u) \right\} $$

where

$$ (3) \quad a(u, v) = \int_{\Omega} (\nabla u \cdot \nabla v + a_o uv) dx $$

$$ (4) \quad b(u) = \int_{\Omega} f u dx + \int_{\Gamma_2} g u dx $$

$$ (5) \quad U(u_o) = \{u : \Omega \to \mathbb{R}, \text{ smooth and such that } u = u_o \text{ on } \Gamma_I \} $$

It is also shown that it makes sense to look for $u$ in the form

$$ (6) \quad u(x) = \sum_{i=1}^{NS} u_i w_i(x), \quad x \in \Omega $$

where $\{w_i\}_{i=1}^{NS}$ span a space $U_h$ and are chosen such that $U_h$ is a good approximation of $U$. In this chapter we consider the simplest approximation. Let $q$ be a triangulation of $\Omega$ i.e.

$$ (7) \quad q = \{ T_i ; i = I, \dots, NT \} \quad U_h = \bigcup_{i=1}^{o} T_i, \quad T_i \cap T_j = \emptyset \quad i \ne j $$

where $T_h$ is a good polygonal approximation of $\Omega$.
For notational convenience we assume that $\Omega = \Omega_h$. Then take

$$ (8) \quad U_h = \{ u \in C^0(\Omega) \mid u|_{T^+} \text{ affine} \} \quad (C^\infty(\Omega) = \text{space of continuous functions}) $$

Note that any $u \in U_h$ is completely determined by its values

$u_i = u(q_i^i)$ on the nodes $\{q_i^i\}_{i=1}^N$ of the triangulation.

Let $\{w_i\}_{i=1}^N$ be the functions of $U_h$ such that

$$ (9) \quad w_i(q_j^i) = 0 \quad i \neq j $$

$$ (10) \quad w_i(q_i^i) = I $$

$\{w_i\}$ is a basis for $U_h$ and any $u \in U_h$ can be written as

$$ (11) \quad u(x) = \sum_i u_i w_i(x) $$

In this frame work problem (2) becomes

$$ (12) \quad \min_{u \in U_h} \{ \frac{1}{2} a(u, u) - b(u) \} $$

Assume that these are NS - MM nodes on $\Gamma_I$ and that they are numbered last, then:

$$ (13) \quad U_h(u_o) = \{u \in U_h | u|_{\Gamma_I} = u_o\} = \left\{ u \mid u = \sum_{i=1}^{MM} u_i w_i + \sum_{i=MM+1}^{NS} u_i w_i \right\} $$

Now it appears that the unknowns in (12) are $\{u_i\}_{i=1}^N$;

therefore (12) is a minimization problem in $R^N$.

## 2. Solution of (12) by the SDM

With the notations of the preliminaries

$$ (14) \quad Z = (u_1, u_2, \dots, u_{MM}) $$

$$ (15) \quad F(Z) = \frac{1}{2} \int_{\Omega} (|\nabla u|^2 + a_0 u^2) dx - \int_{\Omega} f u dx - \int_{\Gamma_2} g u dx $$
and

$$ (16) \quad u = \sum_{i=1}^{M} u_i w_i + \sum_{MM+1}^{N} u_i^T w_i $$

Therefore

$$ (17) \quad \frac{\partial F}{\partial z_i} = \int_{\Omega} (\nabla u \cdot \nabla (\frac{\partial u}{\partial u_i}) + a_o u \frac{\partial u}{\partial u_i}) dx - \int_{\Gamma} g \frac{\partial u}{\partial u_i} dx $$

From (16) : $\frac{\partial u}{\partial u_i} = w_i$ ; hence

$$ (18) \quad \frac{\partial F}{\partial z_i} (Z) = a(u, w_i) - b(w_i) \quad i=1, \dots, M $$

**Exercise 1:** (see Exercise 2, Chapter 1)

Show that $\lambda^i$ in step b of the SDM is given by solving

$$ (19) \quad a(u^i + \lambda d^i, d^i) = b(d^i), \quad d^i(x) = - \sum_{j=1}^{M} \frac{\partial F}{\partial z_j} (u^i) w_j(x) $$

Therefore the SDM applied to (12) is **2.1 Algorithm I**

$$ 0. \quad \underline{\text{set}} \ u^o = \sum_{j=MM+1}^{N} u_j^o w_j; \ \underline{\text{set}} \ i = 0 $$

$$ (20) \quad \underline{\text{compute}} \quad 1. \quad d^i = -a(u^i, w_j) + b(w_j) \quad \hat{\alpha} = 1, \dots, MM $$

$$ \underline{\text{compute}} \quad 2. \quad \alpha^i = a(u^i, d^i) - b(d^i) \quad \text{where} \\ d^i = \sum_{j=1}^{MM} d_j^i w_j $$

$$ \text{and compute} \quad \beta^i = a(d^i, d^i) $$

$$ \text{and set} \quad \lambda^i = -\alpha^i / \beta^i $$

$$ \begin{array}{l} 3. \quad \underline{\text{Set}} \quad u^{i+1} = u^i + \lambda^i d^i \quad (\text{i.e. } u_j^{i+1} = u_j^i + \lambda^i d_j^i, \quad \hat{\alpha} = 1, \dots, MM); \\ \qquad u_j^{i+1} = u_j^o \quad (j = MM+1, \dots, NS), \quad \underline{\text{set}} \ i = i+1 \ \underline{\text{and go back to step 1.}} \end{array} $$
## 2.2. Implementation of algorithm 1

As we shall see all computations will be done on each triangle.
Therefore we need the following information:

- Coordinates of the nodes: arrays $Q(1,J)$, $Q(2,J)$, $J = 1,NS$

- Type of each node: array NG(J) $J = 1, NS$

$NG(J) = 0$ if node J is not on the boundary of $\Omega$

$NG(J) = 1$ belongs to $\Gamma_1$

$NG(J) = 2$ belongs to $\Gamma_2$

- We must be able to know for each triangle what are the numbers of the nodes which are vertices of the triangle:

- array $ME(I,K)$; $I = 1,2,3$; $K = 1, NT$.

As an Example $ME(1,26) = 5$, $ME(2,26) = 7$, $ME(3,26) = 8$ means that triangle no 26 has for vertices the nodes no 5,7,8.

The reader will find in Annex two subroutines SQUARE and CERCLE which triangulate a square and a ring and build the arrays NG, and ME.

### Exercise 2:

Add the center point to the nodes in CERCLE to obtain a triangulation subroutine for the circle. □

<figure><img src="image_6.png" /><figcaption>Figure 1:</figcaption></figure>

Examples of a triangulation with 8 nodes and 8 triangles, numbered.
### Exercise 3:

Construct from SQUARE a triangulation for a Nozzle by moving the nodes of the square, by affine transformations. □

Now we go back to algorithm 1

$$ \text{Computation of } I = \int_{\Omega} \nabla\phi\nabla\psi dx; \quad \phi, \psi \in \mathcal{U}_{\Omega} $$

Since $\Omega = \bigcup_{1}^{NT} T_k$, it suffices to show

how to compute

$$ (24) \quad I_k = \int_{T_k} \nabla\phi\nabla\psi dx $$

Let

$$ (25) \quad \phi_{\partial \Omega} = \phi(q_{\Omega}) $$

then, $\phi$ being linear on $T_k$

$$ \phi(x) = \partial_{\phi_1} x_1 + \partial_{\phi_2} x_2 + \partial_{\phi_3} x_3 \quad \forall x = (x_1, x_2) \in T_k $$

where $\partial_{\phi_1}$, $\partial_{\phi_2}$, $\partial_{\phi_3}$ are constants determined by (25), i.e.

$$ \partial_{\phi_1} q_1^1 + \partial_{\phi_2} q_2^1 + \partial_{\phi_3} q_3^1 = \phi_1 $$

$$ \partial_{\phi_1} q_1^2 + \partial_{\phi_2} q_2^2 + \partial_{\phi_3} q_3^2 = \phi_2 $$

$$ \partial_{\phi_1} q_1^3 + \partial_{\phi_2} q_2^3 + \partial_{\phi_3} q_3^3 = \phi_3 $$
Therefore

$$
(26) \quad \nabla \phi = \left[ \begin{array}{l} \partial \phi_1 = [(\phi_1 - \phi_3)(q_2^2 - q_2^3) - (\phi_2 - \phi_3)(q_1^1 - q_2^3)]/c_2 \\ \partial \phi_2 = -[(\phi_1 - \phi_3)(q_1^2 - q_1^3) - (\phi_2 - \phi_3)(q_1^1 - q_2^3)]/c_2 \end{array} \right.
$$

with

$$
(27) \quad c_2 = (q_1^1 - q_1^3)(q_2^2 - q_2^3) - (q_1^2 - q_1^3)(q_2^1 - q_2^3)
$$

Similarly for $\nabla\psi$. Therefore

$$
(28) \quad I = \sum_{k=1}^{N} I_k \quad \text{with} \quad I_k = \frac{(\partial \phi_1 \partial \psi_1 + \partial \phi_2 \partial \phi_2)}{|c_2|^2}
$$

because $\nabla\phi\nabla\psi$ is constant on $T_k$ and $|c_2|$ is twice the
area of $T_k$

We give a listing of a subroutine, FEM1A, which computes

$$
\mathrm{XA} = \int_{\Omega} (a \nabla \psi \cdot \nabla \phi + a_{0} \psi \phi) dx
$$

the second term being computed by (30) below.

SUBROUTINE FEMIA(A, AO, Q, ME, NS, NT, U, V, XJ)
C COMPUTES XJ=SUM ON OMEGA OF A*GRADU*GRADV*AO*U*V
DIMENSION Q(2,1),ME(3,1),U(1),V(1),Q1(3),Q2(3) , UK(3),VK(3),
1 DU(2),DV(2)
XJ=0.
DO 1 K=1,NT
C PREPARATION ON TRIANGLE K
DO 2 J=1,3
MEKJ=ME(J,K)
Q1(J)=Q(1,MEKJ)
Q2(J)=Q(2,MEKJ)
UK(J)=U(MEKJ)
2 VK(J)=V(MEK)
C COMPUTATION OF GRADU AND GRADV
CALL GRADU(UK,Q1,Q2,DU,C2)
CALL GRADV(VK,Q1,Q2,DV,C2)
IF(C2.LT.O.)C2=-C2
### C FORMULA 2.28 AND 2.30

$$1 XJ=XJ+A*(DU(1)*DV(1)+DU(2)*DV(2)*C2*0.5$$

$$1 * AO* ((UK(1)+UK(2)+UK(3)*(VK(1)+VK(2)+VK(3)))+UK(1)*VK(1)+UK(2)*2VK(2)+UK(3)*VK(3))*C2/24.$$

RETURN

END

$$\text{Computation of } XB = \frac{\int_{\Omega} f\psi dx, \quad \psi \in U_h}{\text{Similarly } XB \text{ is computed from } I_k = \int_{T_k} f\psi dx}$$

If $f$ is constant on $T_k$, then

$$(29) \quad l_k = f(\psi(1) + \psi(2) + \psi(3)) (c_2/2)/3$$

If $f$ is linear on $T_k$ and takes the values $f(1)$, $f(2)$, $f(3)$ on the 3 vertices of $T_k$, then it is easy to show that

$$(30) \quad l_k = \frac{1}{2} c_2 \left[ \begin{aligned} & (f(1) + f(2) + f(3))(\psi(1) + \psi(2) + \psi(3)) + \\ & f(1)\psi(1) + f(2)\psi(2) + f(3)\psi(3) \end{aligned} \right] / 12.$$

If $f$ is not affine $l_k$ may be approximated by a *quadrature formula*

Subroutine FEM1BF will compute $XB$ according to (29) or (30).

SUBROUTINE FEM1BF (Q, ME, NS, NT, V, F, XB, KSW)

COMPUTES SUM ON OMEGA OF F*V BY 2.29 IF KSW=1 OR 2.30 IF KSW=2 RESULT IN XB.

DIMENSION Q(C,1), ME(3,1), V(1), F(1), QK(2,3), VK(3), FK(3)

XB=0.

DO 1 K=1,NT

PREPARATION ON TRIANGLE K

DO 2 J=1,3

MEKJ=ME(J,K)

DO 3 I=1,2

3 QK(1,J)=Q(1,MEKJ)

VK(J)=V(MEKJ)

2 FK(J)=F(MEKJ)

$$C2=(QK(1,1)-QK(1,3))*(QK(2,2)-QK(2,3))-(QK(1,2)-QK(1,3))*(QK(2,1)-QK(2,3))$$

1

IF(C2<0,LT,0.)C2=-C2

IF(KSW.EQ.2)XJ=XJ+((FK(1)+FK(2)+FK(3)*(VK(1)+VK(2)+VK(3))+FK(1)*VK(1)+FK(2)*VK(2)+FK(3)*VK(3))*C2/24.

1

IF (KSW.EQ.1) XJ=XJ+F(K)*(VK(1)+VK(2)+VK(3))*C2/6.

1 CONTINUE

RETURN

END
Computation of $XB = \int_{\Gamma_2} g\psi d\Gamma$, $\psi c U_h$

Similarly one computes $I_e = \int_{C_e} g\psi d\Gamma$ where $C_e$ is the eth straight side of $\Gamma_2$. Let $\psi(1), \psi(2)$ be the values of $\psi$ at the two end points of $C_e$ of coordinates $q^1 = (q_1^1, q_2^1)$ and $q^2 = (q_1^2, q_2^2)$. If $g$ is constant on $C_e$ then

$$ (31) \quad I_e = g(\psi(1) + \psi(2)) \frac{\|q^2 - q^1\|}{2}. $$

If $g$ is linear on $C_e$

$$ (32) \quad I_e = [(g(1) + g(2))(\psi(1) + \psi(2)) + \psi(1)g(1) + \psi(2)g(2)] \\ \qquad \|q^2 - q^1\| / 6. $$
SUBROUTINE PEM18G (Q, ME, NG, NS, NT, V, G, XB, KSV)

C COMPUTES XB = SUM ON GAMMA2 OF G•V ACCORDING TO 2.31 IF KSV=1 OR 2.32

C IF KSV=1 THE VALUES OF G SHOULD BE STORED IN THE ORDER IN WHICH THEY

C AREAR WHEN THE TRIANGLES ARE SCANNED. IF KSV=2 THEN G IS STORED ACCORDING

C TO THE NUMBER OF THE NODE. THEREFORE ONE MAY SAVE COMPUTER TIME BY

C NUMERING THE NODES OF GAMMA2 FIRST.

DIMENSION Q(2,1), ME(3,1), NG(1), V(1), G(1), Q1(3), Q2(3), NGK(3), VK(3)
XB=0
NK=1
DO=1 K=1,NT
NX=0
DO=2 J=1,3
MEK=ME(J,K)
Q1 (J)=Q(1,MEK)
Q2 (J)=Q(2,MEK)
NGX=NG(MEK)
NX=NX+NGX
NGK(J)=NGX

2 VK (J)=V(MEK)

C FINDS IF THE TRIANGLE IS ON GAMMA2

IF (NX.LE.2)GO TO 1

C FINDS THE BOUNDARY NODES J2,J3

DC 3 J=1,3
IF (NGK(J).EQ.0)J1:=J
IF (NGK(J).EQ.2)J2:=J

3 CONTINUE

J3=6-J1-J2

C COMPUTES CL

CL=((Q1(J2)-Q1(J3))^2+(Q2(J2)-Q2(J3))^2)^2(0.5)

C FORMULA 2.31

IF (FSV.EQ.1)XB=XB+G(NK)*(VK(J3)+VK(J2))*CL/2.

C FORMULA 2.32

1F (KSV.EQ.2)XB=XB+((G(J3)+G(J2))*(VK(J3)+VK(J3))+G(J2)*VK(J2)+G(J3)
*VK(J3))*CL/6.

NK=NK+1

1 CONTINUE

RETURN

END
### Exercise 2:

Program algorithm 1 and test it with $\Omega$ = Ring (1,5)
$\phi_0 = 0, a_0 = 0, f = 1, g = 1, \Gamma_1$ = circle [(0,0),1]
$\Gamma_2$ = circle [(0,0), 5] □

### Exercise 3:

1) Same as exercise 2 but with $\Omega = [0,1]^2$, $\Gamma_1 = \partial\Omega$, $\Gamma_2 = \phi$
$\phi_0 = 0, a_0 = 1, f(x_1,x_2) = (\frac{1}{2} + 2\pi^2) \sin\frac{x_1}{2\pi} \sin\frac{x_2}{2\pi}$
Check that the exact solution is: $\phi = 2\pi^2 \sin\frac{x_1}{2\pi} \sin\frac{x_2}{2\pi}$.

2) Take 3 different values of $h$ for the triangulation of $\Omega$
compute $\epsilon_h = \int_{\Omega} |\phi_h - \phi|^2 dx$ and $\epsilon_h^1 = \int_{\Omega} |\nabla(\phi_h - \phi)|^2 dx$;
plot $l_h(\epsilon_h)$, $l_h(\epsilon_h^1)$ versus $l_v(h)$; and find $\alpha$ and $\alpha^1$
such that $\epsilon_h \sim kh^\alpha$ $\epsilon_h \sim k^{-\alpha}h^\alpha$ □

### Exercise 4:

Write the CGM algorithm for problem (12). □

### Exercise 5:

Same as 1) of exercise 2 but with $f(x_1, x_2) = \operatorname{dirac}$ function at $(\cdot, \cdot)$ □
### Remark:

An industrial version of the CCM should not use FEM1A FEM1BF FEM1BG. It takes less time to construct at once $a(u, w_i) - b(w_i)$.

This is done in the subroutine FEM1CG (Finite Element Method P1 with Conjugate Gradient)

a listing of which can be found in appendix. In the sequel we shall use FEM1CG only.

## 3. Solution of (12) by resolution of the linear system

Alternatively (12) is equivalent to

$$ (33) \quad a(u, w_j) = b(w_j) \quad j = 1, \dots, m \text{ mm} $$

$$ (34) \quad a = \sum_{l=1}^{m} u_i w_j + \sum_{mm+1}^{NS} u_{0i} w_i $$

Let

$$ (35) \quad a_{ij} = a(w_i w_j) = \int_{\Omega} (\nabla w_i \cdot \nabla w_j + w_i w_j) dx $$

$$ (36) \quad b_j = b(w_j) = \int_{\Omega} f w_j \, dx + \int_{\Gamma_2} g w_j \, dx $$

then (33) (34) become the linear system (with self explanatory notations)
$$ (37) \quad AU = B $$

A is a positif definite matrix; therefore (37) can be solved by the method of Choleski. To compute A and B see (24) - (32).

This method is faster when many problems (12) with different $f$ and $g$'s must be solved because only $B$ changes and $A$ can be factorized once and for all. However it requires more computer memory. A listing of an efficient subroutine for solving (12) by this method is given in annex. The subroutine's name is

FEMICF (Finite Element method Pl with

Choleski Factorization)

to be used with ASTOF, DPQH, MCHB and a triangulation subroutine. The construction of the subroutine is explained in Chapter 3.

### Exercise 6:

Same as exercise 2 but with FEMICF □

### Remark

There is a simple way (but mathematically not quite correct) to take into account the boundary condition $u_i = u_{0i}$, instead of writing (34) it is easier to write
$$ (38) \quad u = \sum_{i=1}^{NS} u_i w_i $$

and (35) (36) but now in (37) replace

$a_{ii}$ by $10^{20}$ and $b_i$ by $10^{20} a_{0i}$.
RELEASE 2.0

FEMICG

DATE = 77060

15/33/21

SUBROUTINE FEMICG(A, AO, F, G, KSW, KSV, Q, ME, NG, NS, NT, U, KMAX, EPS, ISW, H, 1DJ)

C COMPUTES U SOLUTION OF -D(A*DJ) + AO*U = I IN OMEGA, U ON GAMMA= UO DU/DN ON

C GAMMA 2 = G, BY LAGRANGIAN FINITE ELEMENTS P1. THE SYSTEM IS SOLVED BY
THE CONJUGATE GRADIENT METHOD.

C PARAMETERS. A, AO

F ARRAY (NT) IF KSW=1, (NS) IF KSW=2

G ARRAY (NUMBER OF SIDES OF GAMMA2) IF KSV=1
(NS) IF KSV=2, THEN G IS GIVEN ON THE NODES

Q(3,NS) COORDINATES OF NODES

ME(3,NT) CONTAINS NUMBERS OF NODES OF TRIANGLES

NG(NS)=0 IF INTERIOR NODES. =1 IF DIRICHLET NODE. =2 IF NEUMANN NODE.

NS NUMBER OF NODES

NT NUMBER OF TRIANGLES

KMAX NUMBER MAX OF ITERATIONS OF CONJUGATE GRADIENT.

EPS PRECISION OF CONJUGATE GRADIENT REQUIRED

U MUST BE INITIALIZED TO 0. INSIDE AND UO ON GAMMA 1.

DIMENSION F(1),G(1),Q(2,1),ME(3,1),NG(1),U(1),H(1),DJ(1)

GAMA=0.

DO 1 KK=1,KMAX

CALL GRADJ(Q,ME,NS,NT,NG,A, AO, F, G, KSW, KSV, DJ, U)

DJ2=0.

XJ=0.

GAM=0.

DO 4 I=1,NS

4 GAM=GAM+DJ(I)**2

IF (KK,GT,1)GAMA=GAM/GAMA

IF (ISW,NE,2)GO TO 5

CALL AHH(Q,ME,NS,NT,NG,A, AO, U, XJ)

XJ=-XJ/2.

5 DO 2 I=1,NS

XJ=XJ+DJ(I)*U(I)

H(I)=-DJ(I)*GAMA*H(I)

2 DJ2=DJ2-DJ(I)*H(I)

CALL AHH(Q,ME,NS,NT,NG,A, AO, H,HA2)

XLAM= DJ2/HA2

DO 3 I=1,NS

3 U(I)=U(I)+H(I)*XLAM

IF (ISW,EQ,2)PRINT 100, XJ, DJ2, XLAM, GAMA

100 FORMAT (IX,') COUT = ',E13.6,' DJ**2 = ',E9.2,' LAM = ',E9.2,

1 GAMA = ',E9.2)

IF (DJ2,LE,EPS)RETURN

GAMA=GAM

1 CONTINUE

RETURN

END

RELEASE 2.0

FEMICG

DATE = 77060

15/33/21

EFFECT* NOTERM, ID, EBCDIC, SOURCE, NOLIST, NODECK, LOAD, NOMAP, NOTEST

EFFECT* NAME = FEMICG , LINECNT = 50

* SOURCE STATEMENTS = 29, PROGRAM SIZE = 1510

* NO DIAGNOSTICS GENERATED
RELEASE 2.0

GRADJ

DATE = 77060

15/33/21

SUBROUTINE GRADJ(Q,ME,NS,NT,NG,A,AO,F,G,KSW,KSV,DJ,U)

C COMPUTES DJ/DWI WHERE J= SUM/OMEGA OF (A*DU*DU/2*AO*U*U/2 -F*U)

C *SUM/GAMMA2 OF G*U. WHERE U= SIGMA(UI*WI)
DIMENSION Q(2,1),MF(3,1),F(1),G(1),NG(1),DJ(1),U(1),DU(2),MEK(3),
1 Q1(3),Q2(3),NGK(3),UK(3),DW(2),VK(3),WJ(3)

C PREPARATION

L1=1
L2=2
X2=.5
X3=1./3.
X6=X2*X3
X12=X6*X2
DO 3 I=1,NS

3 DJ(I)=0.
KF=0

DO 1 K=1,NT

C PREPARATION ON TRIANGLE K

NK=0
DO 2 J=1,3
MEKJ=ME(J,K)
MEK(J)=MEKJ
Q1(J)=Q(1,MEKJ)
Q2(J)=Q(2,MEKJ)
NGK(J)=NG(MEKJ)
UK(J)=(I)(MEKJ)
VK(J)=UK(J)*AO
IF(KSW.EQ.2)VK(J)=VK(J)-F(MEKJ)

? NK=NK*NGK(J)

C FINDS THE INTERIOR NODE IF K IS A BOUNDARY TRIANGLE

KSVV=0
IF(NK-3)11,13,13

13 KF=KF+1
KSVV=KSV
DO 14 J=1,3
J1=J

IF(NGK(J).LE.0)GO TO 11

14 CONTINUE

11 CONTINUE

C COMPUTES SUM(A*DU*DWI)

CALL GRADU(UK,Q1,Q2,DU,AIR,L1,DN)
DO 1 J=1,3
IF(NGK(J).EQ.1)GO TO 1
WJ(1)=(J-2.)*(J-3.)*X2
WJ(2)=(J-1.)*(J-3-J)
WJ(3)=(J-1.)*(J-2.)*X2
MEKJ=MEK(J)
CALL GRADU(WJ,Q1,Q2,DW,AIR,L2,DN)
DJK=A*(DU(1)*DW(1)+DU(2)*DW(2))

C COMPUTES SUM(U*WI-F*WI)

IF(KSW-1)6,5,6

5 DJK=DJK-F(K)*X3

6 DO 10 I=1,3

10 DJK=DJK+(1.*WJ(I))*VK(I)*X12

C COMPUTES AREA OF TRIANGLE

IF(AIR.LT.0.)AIR=-AIR
DJK=DJK*AIR

C COMPUTES SUM G*WI ON GAMMA2 IF NECESSARY

4 IF(NGK(J).NE.2) GO TO 16

IF(KSVV-1)16,8,9

8 J2=3-J1-J
XL=((Q1(J)-Q1(J2))**2+(Q2(J)-Q2(J2))**2)**X2
DJK=G(KF)+XL*X2+DJK
GO TO 16

9 DJK=G(MEKJ)*XL*X3+G(MEK(J2))*XL*X6+DJK

16 DJ(MEKJ)=DJ(MEKJ)+DJK

1 CONTINUE
IF(KSW.NF.0)RETURN
DO 15 I=1,NS

15 DJ(I)=DJ(I)+F(I)
RETURN
END
1 RELEASE 2.0

AHH

DATE = 77060

15/33/21

SUBROUTINE AHH(Q,MF,NS,NT,NG,A,A0,H,HA2)

C CALCUL SUM/OMEGA (A*DH*DH+A0*H*H)

DIMENSION Q(2,1),MF(3,1),NG(1),H(1),Q1(3),Q2(3),HK(3),DHK(2)

L1=1

A12=A0/12.

X2=.5

HA2=0.

DO 1 K=1,NT

DO 2 J=1,3

MEK=ME(J,K)

Q1(J)=Q(1,MEK)

Q2(J)=Q(2,MEK)

2 HK(J)=H(MEK)

CALL GRADU(HK,Q1,Q2,DHK,AIR,L1,DN)

IF(AIR,LT,0.)AIR=-AIR

1 HA2=HA2*((DHK(1)**2+DHK(2)**2)*A*((HK(1)•HK(2)•HK(3)))**2•HK(1)**2

1 •HK(2)**2•HK(3)**2)*A12)*AIR

RETURN

END

1 RELEASE 2.0

GRADU

DATE = 77060

15/33/21

SUBROUTINE GRADU(U,Q1,Q2,DU,AIR,JSW,DN)

C COMPUTES GRADIENT U ON A TRIANGLE OF COORDINATES OF VERTICES Q1,Q2.

C IF JSW NOT=1 AIR IS INPUT.

DIMENSION U(3),Q1(3),Q2(3),DU(2)

Q113=Q1(1)-Q1(3)

Q213=Q2(1)-Q2(3)

Q123=Q1(2)-Q1(3)

Q223=Q2(2)-Q2(3)

U13=U(1)-U(3)

U23=U(2)-U(3)

IF (JSW=1) 1,2,1

? AIR=J.5*(Q113*Q223-Q123*Q213)

DN=.5/AIR

1 CONTINUE

DU(1)=(U13*Q223-U23*Q213)*DN

DU(2)=(Q113*U23-Q123*U13)*DN

RETURN

END
G1 RELEASE 2.0

SQUARE

DATE = 77060

15/33/21

SUBROUTINE SQUARE(M,N,NS,Q,NT,MF,NG,MUD,NF,JSW)

C TPIangsULATeS A SoUARe By 3 FAmIlIES OF eQUIDISTANT LINES

C M = NUMBER OF NODES ON THE HORIZONTAL SIDE

C N = VERTICAL

C NS = NUMBER OF NODES

C Q(2,NS) = COORDINATES OF THE NODES

C NT = NUMBER OF TRIANGLES

C MF(3,NT) = NUMBERS OF THE NODES VERTICES OF THE TRIANGLE

C NG(NS) = 0 IF THE NODE IS INTERIOR

C 1 IF THE NODE IS ON THE BOUDARY

C MUD = BANDWIDTH OF THE CHOLESKI MATRIX ( OUTPUT)

C NF = NUMBER OF BOUDARY NODES

C IF JSW = 1 ME,Q,NG ARE PRINTED
DIMENSION Q(2,1),MF(3,1),NG(1)
NS=N*M
NT=2*(N-1)*(M-1)
DO 1 I=1,M
S=I-1.
S=S/(M-1.)
DO 1 J=1,N
R=J-1.
R=R/(N-1.)
K=(J-1)*M*I
NG(K)=0
Q(1,K)=R
1 Q(2,K)=S
DO 2 J=1,M
NG(J)=1
2 NG(M*(N-1)+J)=1
DO 22 J=1,N
NG(N*(J-1)+1)=1
22 NG(N*(J-1)*N)=1
NT2=NT/2
DO 12 K=1,NT2
J=(K-0.1)/(M-1)
I=K-J*(M-1)
MF(1,2*K-1)=I+J*M
ME(2,2*K-1)=I+1+J*M
MF(3,2*K-1)=I+1+J*M+M
L=2*K
MF(1,L)=I+J*M
MF(2,L)=I+1+J*M+M
MF(3,L)=I+J*M+M
12 CONTINUE
MUD=M+2
NF=2*N+2*M-4
IF (JSW,NF,1) RETURN
PRINT 102,NS,NT,MUD,NF

G1 RELEASE 2.0

SQUARE

DATE = 77060

15/33/21

102 FORMAT(1X,' NS = ',.13,' NT = ',.13!) MUD = ',.13,' NF = ',.13)
PRINT 103

103 FORMAT(1X,'FIRST ARRAY IS ME, SECOND IS Q, THIRD IS NG')
PRINT 100, (ME(1,K),ME(2,K),ME(3,K),K=1,NT)

100 FORMAT(10(2X,3I3))
PRINT 101, (Q(1,K),Q(2,K),K=1,NS)

101 FORMAT(6(2X,2E9.2))
PRINT 100 , (NG(K),K=1,NS)
RETURN
END
G1 RELEASE 2.0

CERCLE

DATE = 77060

15/33/21

SUBROUTINE CERCLE (N,M,RO1,RO2,Q,ME,NG,NT,NS,NF,MUD,JW)
C TRIANGULATES A RING BY INTERSECTION OF RAYS AND CERCLES DIVIDED INTO 2
C FROM BOTTOM LEFT TO TOP RIGHT FOR THE TRIANGLES ON THE RIGHT AND
C BY CONTINUITY FOR THE OTHERS.
C INPUT N NUMBER OF NODES ON EACH CIRCLE
C M NUMBER OF CERCLES
C RO1 INNER BOUDARY RADIUS
C RO2 OUTER BOUDARY RADIUS
C JW =2 Q,ME,NG ARE PRINTED OTHERWISE THEY ARE NOT.
CC OUTPUT Q(2,NS) COORDINATES OF THE NODES
C ME(3,NT) NUMBERS OF THE VERTICES OF EACH TRIANGLES
C NG(NS)=0 IF THE NODE IS INTERIOR
C =1 IF THE NODE IS ON THE BOUDARY
C NT NUMBER OF TRIANGLES
C NS NUMBER OF NODES
C NF NUMBER OF BOUDARY NODES
C MUD BANDWIDTH OF THE MATRIX OF CHOLESKI.
DIMENSION Q(2,1),MF(3,1),NG(1)
PIN=2*3.141592653589/N
NS=N*M
NT=2*N*(M-1)
NF=N*2
MUD=N*2
XK=(RO2-RO1)/(M-1)
DO 1 I=1,M
DO 1 J=1,N
NIJ=N*(I-1)+J
Q(2,NIJ)=(RO1+ XK* (I-1))* SIN((J-1)*PIN)
Q(1,NIJ)=(RO1+ XK* (I-1))* COS((J-1)*PIN)
NG(NIJ)=0
K=2*N*(I-1)+2*(J-1)+1
IF (I.EQ.M)NG(NIJ)=1
IF (I.EQ.1)NG(NIJ)=1
IF (K.GT.NT)GO TO 1
ME(1,K)=NIJ
ME(2,K)=NIJ+N
ME(3,K)=NIJ+N+1
ME(1,K+1)=NIJ
ME(2,K+1)=NIJ+N+1
ME(3,K+1)=NIJ+1
IF (J.NE.N)GO TO 1
ME(2,K+1)=NIJ+1
ME(3,K)=NIJ+1
ME(3,K+1)=NIJ-N+1
1 CONTINUE
IF (JW.NE.1)RETURN
PRINT 200,((ME(I,K),I=1,3),K=1,NT)
200 FORMAT(8(2X,3(1X,13)))

1 RELEASE 2.0

CERCLE

DATE = 77060

15/33/21

PRINT201,(I,0(1,I),Q(2,I),NG(I),I=1,NS)
201 FORMAT(4(1X,13,E9.2,2X,E9.2,2X,11))
RETURN
END
# Chapter 3

## Basic Concepts About the Implementation of the Finite Element Method

### Introduction

Our purpose is double. First to introduce the reader to the concepts necessary for the programming of a simple problem by elementary finite elements.

Then to introduce some more sophisticated methods necessary for the solution of large problems (that is problems requiring very much computing time and memory space).

We shall consider only linear problems, and *direct methods* for the solution of the final linear system. The case of non linear problems will be considered later, and the use of iterative methods for solving linear problems was considered in Chapter 2.

Let us describe the main sections

A - Basic concepts.

1. Form of the data necessary for the representation of a triangulation.

2. Stiffness submatrix of an element. Stiffness matrix. Algorithm for the "assembling" of the stiffness matrix.

3. Computation of a stiffness submatrix.

4. Solution of the linear system by the Choleski method - choice of a numbering for the nodes.

B. Computation of the stiffness submatrix for more sophisticated elements.

C. Various methods to compute automatically a triangulation of a domain.
D - Description of two contradictory methods for the solution of large problems.

- The frontal method

- The sub-structure method

A - Basic Concepts

1. Form of the datas necessary for the representation of a triangulation.

Let us consider a simple triangulation $T$.

The following informations are needed:

- The coordinates of the vertices.

- The knowledge of the vertices which are really joined.

<figure><img src="image_7.png" /><figcaption>Fig. 1 Triangulation $T$.</figcaption></figure>

For these purposes, it is convenient to introduce the following array: $M(3 \times n_T)$ where $n_T$ is the number of triangles. We have $n_T = 11$ in the case of Fig. 1, and

$$ M = \begin{pmatrix} 1 & 2 & 6 & 5 & 5 & 3 & 5 & 5 & 4 & 1 & 10 \\ 2 & 3 & 9 & 6 & 6 & 6 & 8 & 8 & 8 & 4 & 5 \\ 5 & 6 & 7 & 9 & 2 & 7 & 9 & 10 & 10 & 10 & 1 \end{pmatrix} = (m_{ij}) \quad 1 \le i \le 3 \\ 1 \le j \le 10 $$

$m_{ij}$ is the number of the $i$th vertex of the $j$th triangle (see Fig. 1).

Let $n_v$ be the number of vertices of the triangulation, $m_{ij}$ is then an integer such that $1 \le m_{ij} \le n_v$.

We introduce then the array $Z = (z_{ij})_{\begin{smallmatrix} 1 \le i \le 2 \\ 1 \le j \le n_v \end{smallmatrix}}$
such that $z_{ij}$ is the i-th coordinate of the j-th vertex.

With these two datas M and Z it is very easy to find
the coordinates of the vertices of some triangle, and then to
compute for example the area of this triangle.

As additional information about the triangulation, one
may require to know if a vertex belongs to the boundary which
can be specified by the data of a simple vector

$$L = (\ell_i)_{1 \le i \le n_v}$$

where $\ell_i = 1$ if the vertex i belongs to the boundary of the domain, and $\ell_i = 0$ otherwise.

In the following, we shall assume that the arrays M, Z
and eventually L are given.

In fact the knowledge of the quantities M and Z is sufficient
to draw the triangulation T, but it is sometimes easier to know
the vector L.

### Exercise 1:

Draw the triangulation T' corresponding to the datas

$$M = \begin{pmatrix} 2 & 2 & 6 & 3 & 4 & 4 & 1 & 5 \\ 7 & 5 & 5 & 6 & 6 & 9 & 7 & 7 \\ 5 & 8 & 8 & 8 & 5 & 5 & 9 & 9 \end{pmatrix}$$

$$Z = \begin{pmatrix} 0. & 1. & 1. & 0. & 0.5 & 0.5 & 0.5 & 1. & 0. \\ 0. & 0. & 1. & 1. & 0.5 & 1. & 0. & 0. & 0.5 \end{pmatrix}$$

### Exercise 2:

Consider the domain $\Omega = ]0,1[x]0,1[$, and the uniform mesh
obtained by dividing $\Omega$ in N x N equal squares. Write the FORTRAN
program which computes the datas M and Z for the triangulation
obtained by dividing each square in two triangles by one of its diagonals.

(Note that in M, the order in which the 3 vertices of a triangle appear has no importance. However it can be convenient to order them anti-clockwise)

## Other Types of Triangulations and Finite Elements

Let us consider the case of a "triangulation" made with quadrilateral elements (Fig. 2).

<figure><img src="image_8.png" /><figcaption>Fig. 2 A Triangulation composed with quadrilateral</figcaption></figure>

In such a case the triangulation will be conveniently represented by the array

$$M = \begin{pmatrix} 1 & 2 & 3 & 5 & 9 \\ 2 & 3 & 4 & 6 & 10 \\ 9 & 10 & 5 & 7 & 7 \\ 8 & 9 & 10 & 10 & 8 \end{pmatrix}$$

(the array Z will be the one you think). Note that the order in the columns of M is now essential, because in a quadrilateral one vertex is only connected to 2 other vertices.

For higher order elements it is convenient to number all the nodes corresponding to degrees of freedom, and not only vertices. As an example for Lagrange finite elements of order 2 (Fig. 3), the array M will be 6 x r, because there are 6 nodes per elements.
$$ M = \begin{pmatrix} 1 & 3 & 3 & 3 & 5 \\ 3 & 11 & 9 & 5 & 7 \\ 13 & 13 & 11 & 9 & 9 \\ 14 & 12 & 10 & 18 & 8 \\ 15 & 14 & 16 & 17 & 18 \\ 2 & 16 & 17 & 4 & 6 \end{pmatrix} $$

<figure><img src="image_9.png" /><figcaption>Fig.3 An example of triangulation with higher degree elements</figcaption></figure>

(Note that the order in the columns of M is always the same, in some sense.)

## 2. Stiffness Submatrix and Stiffness Matrix

Let us assume, for the sake of simplicity, that we want to solve the following variational problem (where $\Omega \subset \mathbb{R}^2$ is a polygonal domain).

$$ (1) \quad \int_{\Omega} (uv + \operatorname{grad} u \operatorname{grad} v) dx = \int_{\Omega} fvdx, \quad \forall v \in H^1(\Omega). $$

Note that (1) is the variational formulation of the Neumann problem introduced by LIONS.

Let $T_h$ be a triangulation of $\Omega$, we introduce the finite element space $V_h$ of the functions which are continuous on $\bar{\Omega}$, the restriction of which to any triangle of $T_h$ is affine. (We shall consider other elements further).

The approximate problem is then to find $u_h \in V_h$ such that

$$ (2) \quad \int_{\Omega} (u_h v_h + \operatorname{grad} u_h \operatorname{grad} v_h) dx = \int_{\Omega} fv_h dx \quad \forall v_h \in V_h. $$

Let $(W_i(x))_{1 \le i \le N}$ ($N$ depending on $h$) be a basis of $V_h$ and
$$ (3) \quad u_h(x) = \sum_{i=1}^{N} y_i W_i(x). $$

We know that the components $(y_i)_{1 \le i \le N} = y \in \mathbb{R}^N$ are solutions of a linear system

$$ (4) \quad Ay = b $$

where $A = (A_{ij})_{1 \le i,j \le N}$ is the so called "stiffness matrix", such that

$$ (5) \quad A_{ij} = \int_{\Omega} (W_i W_j + \operatorname{grad} W_i \operatorname{grad} W_j) dx, \quad 1 \le i, j \le N. $$

The vector $b = (b_i)_{1 \le i \le N} \in \mathbb{R}^N$ is such that

$$ (6) \quad b_i = \int_{\Omega} f w_i dx, \quad 1 \le i \le N. $$

A simple basis of the space $V_h$ is the following. For $i \in [1,N]$,
$W_i$ is this function of $V_h$, the value of which is 1 on the vertex
number $i$, and 0 on any other vertex of the triangulation. (A
function of $V_h$ is uniquely determined by its values at the vertices
of the triangulation).

However, the use of the basis function $W_i$ and the formula (5) is a very unefficient way to compute the stiffness matrix $A$, because one has to compute several times a quantity relative to one triangle!

To avoid this, we introduce the notion of *stiffness submatrix*.

## Stiffness Submatrix

We remark first that
$$ (7) \int_{\Omega} (u_h v_h + \operatorname{grad} u_h \operatorname{grad} v_h) dx = \sum_{K} \int_{K} (u_h v_h + \operatorname{grad} u_h \operatorname{grad} v_h) dx. $$

Let us still denote by K the number of the triangle K;
the numbers of its 3 vertices are $m_{1K}$, $m_{2K}$, $m_{3K}$.

Let $y_i$ (resp. $z_i$) denote the value of $u_h$ (resp. $v_h$) at the vertices of the triangulation, and

$$ \lambda_{\alpha}^{K}(x) \qquad \alpha = 1, 2, 3 $$

the barycentric coordinates associated to the nodes

$$ m_{\alpha K} \qquad \alpha = 1, 2, 3. $$

The value $u_h(x)$ (resp. $v_h(x)$) for $x \in K$ is given by

$$ u_h(x) = \sum_{\alpha=1}^{3} y_{m_{\alpha K}} \lambda_{\alpha}^{K}(x) $$

$$ (\operatorname{resp. } v_h(x) = \sum_{\alpha=1}^{3} z_{m_{\alpha K}} \lambda_{\alpha}^{K}(x)). $$

We define then the stiffness submatrix $A^K$ of K as the 3 x 3
matrix defined by

$$ (8) \quad A_{\alpha\beta}^{K} = \int_{K} (\lambda_{\alpha}^{K} \lambda_{\beta}^{K} + \operatorname{grad} \lambda_{\alpha}^{K} \operatorname{grad} \lambda_{\beta}^{K}) dx, \quad 1 \le \alpha, \beta \le 3. $$

From this definition follows the identity

$$ (9) \quad \int_K (u_h v_h + \operatorname{grad} u_h \operatorname{grad} v_h) dx = \sum_{\alpha, \beta=1}^{3} y_{m_{\alpha K}} z_{m_{\beta K}} A_{\alpha\beta}^{K}. $$

Then, from the definition (5) of A, it follows that
$$ (10) \quad \sum_{i,j=1}^{N} y_i z_i A_{ij} = \sum_{K \in T_h} \sum_{\alpha, \beta=1}^{N} y_{m_{\alpha K}} z_{m_{\beta K}} A_{\alpha \beta}^{K} $$

This formula being true for any y, z ∈ R<sup>N</sup>, we see that the stiffness matrix A can be computed from the submatrices A<sub>K</sub><sup>K</sup>, for K ∈ T<sub>h</sub>.

We can check that the following algorithm achieves the computation of A from the A<sub>K</sub><sup>K</sup>'s (we use a language adapted from programming).

$$ \begin{array}{ll} \text{• \underline{Make} } A = O & (A \text{ is a } N \times N \text{ matrix}). \\ \text{• \underline{For each } } K \in T_h & \\ \quad \text{• \underline{compute }} A^K & \\ \text{(11) } \quad \text{• \underline{For }} \alpha, \beta = 1, 2, 3, \underline{make} & \\ \quad \quad A_{m_{\alpha K}, m_{\beta K}} = A_{m_{\alpha K}, m_{\beta K}} + A_{\alpha \beta}^K & \end{array} $$

The algorithm (11) is called "assembling" of the matrix A; its advantage is that you only need to compute the quantities relative to one element once.

In Fig. 4, we represent the evolution of the matrix A during the assembling. The stiffness submatrices of the elements a, b, c, d, will be denoted by $\begin{pmatrix} a & a & a \\ a & a & a \\ a & a & a \end{pmatrix}, \begin{pmatrix} b & b & b \\ b & b & b \\ b & b & b \end{pmatrix}, \dots$

because we are not interested by the actual values of these coefficients.
<figure><img src="image_10.png" /><figcaption>Fig. 4 An example of assembling</figcaption></figure>

### Exercise 3:

Write the algorithm (11) in FORTRAN. (The arrays A(N, N), AK(3, 3), M(3, NT) are supposed to have been already declared by an order DIMENSION).

### Remark 1:

Let us consider now the case of a Dirichlet boundary condition:

$$ (12) \quad u_h = 0 \quad \text{on } \Gamma $$

In the variational foundation (2), $v_h$ will satisfy the same boundary condition. Then the new stiffness matrix (the one to use in the linear system (4)) will be obtained from the one given by the assembling (11) by suppression of all the lines and columns relative to boundary nodes.

Of course, in a practical way, it is better to adapt the algorithm (11) to this new situation. For this purpose we introduce a vector noted iRA (in Fortran) such that for $1 \le j \le n$:
iRA(j) = 0 if the vertex j belongs to the boundary
iRA(j) is a new number (ε[1, n_w] where n_w is the number of
internal vertices) otherwise.

The new assembling algorithm will be the following:

• Make A = 0

• For each K ∈ T_h:

- compute A^K

- For α, β = 1, 2, 3 compute

n_α = iRA(m_α^K)

n_β = iRA(m_β^K)

(13)

if n_α = 0 or n_β = 0, then do nothing;

otherwise make

A^{n_{\alpha}}, n_{\beta} = A^{n_{\alpha}}, n_{\beta} + A^{K}_{\alpha\beta}.

**Exercise 4** (analogous to exercise 3):

Write the algorithm (13) in FORTRAN.

**Remark 2:**

Note that the array iRA can be conveniently computed
from the data L defined in §1 by the following FORTRAN
instructions

j = 0
DØ 1 i=1, NV
iRA(i)=0
iF(L(i)) 2, 2, 1
2 j=j+1
iRA(i)=j
1CØNTINUE .
**Remark 3:**

Note that all these assembling algorithms can be easily extended to the case of other finite elements. The stiffness submatrix will only have a different size. As an example, the dimensions of $A^K$ for a Lagrange finite element of type 2, will be 6 by 6 (instead of 3 x 3).

### 3. Computation of a Stiffness Submatrix

Let i, j, k denote the number of the vertices of a triangle K. We shall consider a reference element $\hat{K}$, the vertices of which are (0,0), (1,0), (0,1). We consider the affine mapping F which associates i to (0,0), j to (1,0), and k to (0,1)(Fig. 5).

<figure><img src="image_11.png" /><figcaption>Fig. 5 The Affine Mapping F</figcaption></figure>

Let $(x_i, y_i)$, $(x_j, y_j)$, $(x_k, y_k)$ denote the coordinates of the vertices i, j, k. We can easily check that

$$F(\xi, \eta) = B(\xi) + c$$

where B is 2 x 2 matrix
$$
(14) \quad B = \begin{pmatrix} x_j - x_i & x_k - x_i \\ y_j - y_i & y_k - y_i \end{pmatrix}
$$

and

$$
c = \begin{pmatrix} x_i \\ y_i \end{pmatrix} .
$$

Infact ξ, η are connected to the barycentric coordinates
in the following way:

$$
(15) \quad \lambda_i^K = 1 - \xi - \eta; \quad \lambda_j^K = \xi; \quad \lambda_k^K = \eta .
$$

The matrix B is called Jacobian matrix, and we have

$$
(16) \quad \left| \det B \right| = \frac{\operatorname{mes}(K)}{\operatorname{mes}(\hat{K})} = 2 \operatorname{mes}(K)
$$

where mes(K) denotes the area of the triangle K.

Moreover, if the nodes i, j, k are ordered anti clockwise,
det B > 0, and the absolute value is not necessary in (16).

Let now v(x, y) be a function defined on K, and denote
by v(ξ, η) the corresponding function defined on K, such that

$$
\hat{v}(\xi, \eta) = v(F(\xi, \eta)).
$$

One can easily check that

$$
\left. \nabla v \right|_{\xi, \eta} = B^T \left. \nabla v \right|_{F(\xi, \eta)}
$$

Then the following formula holds:

$$
(17) \quad \int_K \nabla u \cdot \nabla v \, dxdy = \int_{\hat{K}} ((B^T)^{-1} \nabla u) \cdot ((B^T)^{-1} \nabla v) \cdot |\det B| \, d\xi \, d\eta
$$

This formula can be used for the computation of the term

$$
(18) \quad \int_K \nabla \lambda_m^K \cdot \nabla \lambda_{m\beta}^K \, dxdy
$$
which occurs in the definition (8) of the stiffness submatrix
$A_i^K$. In fact the formula (15) should be more rigorously written

$$
\hat{\lambda}_i^K = 1 - \xi - \eta; \quad \hat{\lambda}_j^K = \xi; \quad \hat{\lambda}_k^K = \eta.
$$

An easy computation gives

$$
(19) \quad \nabla \hat{\lambda}_i^K = \begin{pmatrix} -1 \\ -1 \end{pmatrix}; \quad \nabla \hat{\lambda}_j^K = \begin{pmatrix} 1 \\ 0 \end{pmatrix}; \quad \nabla \hat{\lambda}_k^K = \begin{pmatrix} 0 \\ 1 \end{pmatrix}.
$$

Then the term to be computed in (18) is completely determined
by (17) and (19). For example one has

$$
\int_K \nabla \lambda_i^K \nabla \lambda_j^K \, dxdy = \frac{1}{2} \left( (B^T)^{-1} \begin{pmatrix} -1 \\ -1 \end{pmatrix} \right) \cdot \left( (B^T)^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} \right) |\det B|
$$

The computation of the term $\int_K \lambda_i^K \lambda_j^K \, dxdy$ which occurs in (8)
is in fact simpler from the formula:

$$
(20) \quad \int_K (\lambda_i!)^\alpha (\lambda_j!)^\beta (\lambda_k!)^\gamma \, dxdy = 2 \operatorname{mes}(K) \frac{\alpha! \beta! \gamma!}{(\alpha+\beta+\gamma+2)!}
$$

which gives immediately the result.

Remark 4:

If the variational formulation (1) changes, then one
has to change in a corresponding way the definition (8) of the
stiffness matrix $A$.

Exercise 5:

Consider the case of a Dirichlet problem: find
$u \in H_0^1(\Omega)$ such that
$$ (21) \quad \int_{\Omega} \operatorname{grad} u \cdot \operatorname{grad} v \, dx = \int_{\Omega} f v \, dx \quad \forall v \in H_0^1(\Omega) $$

instead of the Neumann problem (1). We still use order 1 elements. Compute the stiffness matrix $A^1$ corresponding to the triangle $T_1$, the vertices of which are

$$ (a, b), (a + h, b), (a, b + h). $$

Compute then $A^2$ corresponding to the triangle $T_2$ where $h$ is changed in $-h$.

Is the result depending on $a$ and $b$?

### Exercise 6:

Use the results of the previous exercise and the assembling algorithm to find the 5th line of the stiffness matrix associated to the triangulation given in Fig.6.

<figure><img src="image_12.png" /><figcaption>Fig. 6</figcaption></figure>

(The result proves that this simple finite element method with such a uniform mesh is equivalent to a classical finite difference scheme).
## 4. Solution of the Final Linear System

The computation of the right hand side b occurring in the linear system (4) can be made during the assembling algorithm ((11) or (13)) in the same way that the stiffness matrix A is computed (you can check it as a very useful exercise).

For solving the linear system (4), one can think of using Cholesky's method, because the matrix A is symmetric and positive definite, and in fact it is possible to use it.

However the matrix A is sparse and considerable computing time can be saved by avoiding unnecessary operations on elements which are equal to zero.

The process is very easy when the matrix has a "band-structure" (Cf. Fig.7): For some number w called "band width" one has

$$|i-j| > W \Rightarrow A_{ij} = 0.$$

<figure><img src="image_13.png" /><figcaption>Fig. 7 Band-Structure</figcaption></figure>

For such matrices it is very easy to avoid unnecessary operations, and there exists even a standard IBM scientific program adapted to this case.
The number of operations for a full $N \times N$ matrix (or considered as full) is asymptotically $\frac{N^3}{3}$ operations, while if the matrix has a band width $w$ this number reduces to $\frac{Nw^2}{2}$.

We shall see that the stiffness matrix $A$ has a band structure if the nodes are suitably numbered.

In fact the stiffness matrix $A$ has the following fundamental property:

If $A_{ij} \neq 0$, then there exists an element $K \in T_h$ such that the nodes $i$ and $j$ belong to $K$. In such a case we shall say that the nodes $i$ and $j$ are **neighbours**.

If we want to give to the matrix $A$ a band structure, then it is sufficient that any two neighbours of the trian-gulation have numbers not too different. The worst case is when the first and the last node are neighbours.

<figure><img src="image_14.png" /><figcaption>Fig. 8 The natural way of numbering the nodes for a uniform mesh on a square (here $N = 36$; $w = 6$)</figcaption></figure>
For a uniform mesh on a square (Fig. 8) there is a way of numbering which produces a very small bandwidth: for a n x n square mesh, the bandwidth *w* can be equal to *n*. This gives a number of operations which is $\frac{n^4}{2}$ instead of $\frac{n^6}{3}$. For more complicated geometries the optimal numbering for the nodes looks very much the one of Fig. 9:

<figure><img src="image_15.png" /><figcaption>Fig. 9. - Numbering of the nodes in the case of a general geometry</figcaption></figure>

Note that to evaluate $W$, one has just to compute the maximum of $|i-j|$ for *i* and *j* neighbours. In the case of Fig. 8, the maximum is reached for $i = 10$ and $j = 16$, then $W = 6$.

A simple (and not too bad) way of numbering the nodes is the Cuthill and McKee algorithm:

- give the number 1 to a boundary node

- give the following numbers to its neighbours (2, 3, 4 in the case of Fig. 9)

- give the following numbers to the neighbours of the neighbours (5, 6, 7 in the case of Fig. 9)

- etc. ...
We should now precise that the Cholesky algorithm conserves the band structure of the matrix A. That is, when one computes $A = LL^T$, where L is lower triangular, it is easy to see that L has the same bandwidth W than A.

Note that the application of Gaussian elimination would have the same effect (this because the 2 methods are equivalent in some sense).

### Exercise 7:

Use the induction formula for the computation of the factorized L in Cholesky method to prove that L has the same bandwidth than A.

### Conclusion for part A

Now the reader should be able to program without difficulty a simple finite element program such as the solution of a Dirichlet problem on a circle, approximated by the simplest finite elements. (We advise the reader to use the standard IBM program for the Cholesky band solution). (*)

In the following we shall give some complements. First for more sophisticated elements, then for more complicated geometries, and, at last for large problems.
B - Computation of the Stiffness Submatrix for more Sophisticated Elements

The method which consists (Cf. §A-3) of computing the integrals on a reference element $\hat{K}$ is used most of the times.

When the elements are "straight", that is, when the mapping $\hat{F}: \hat{K} \to K$ is affine, all the integrals (for the computation of the stiffness submatrix $A_K$) can be computed exactly by evaluating them on the reference triangle $\hat{K}$, and by using the formula analogous to (20):

$$ (22) \qquad \int_{\hat{K}} \xi^{\alpha} \eta^{\beta} d\xi d\eta = \frac{\alpha! \beta!}{(\alpha + \beta + 2)!} $$

Of course if the elements are parallelograms the reference element is a square, and the formula is still simple.

However in the more general case of isoparametric elements this is no longer true. One has to use a *quadrature scheme* on the reference element, that is, replace any integral of the type

$$ \int_{\hat{K}} \psi(\hat{x}) d\hat{x} $$

by

$$ \begin{cases} 1. & \psi(\hat{b}_e) \hat{w}_e \\ \lambda & e=1 \end{cases} $$

where the points $\hat{b}_e$ are given in $\hat{K}$, and the coefficients $\hat{w}_e$ are suitably chosen (see LIONS [1]). Such a quadrature scheme can be necessary in any case to evaluate the right hand side $b$. In the following, we shall consider only the case of Lagrange finite elements.

Let us denote by $(p_i(\xi, \eta))_{1 \le i \le n_D}$ the basis functions of the reference finite element $\hat{K}$. Let $v$ be a function of the
approximate space $V_h$ and $(v_i)_{i \le j \le n_v}$ denote the values of $v$ at the $n_v$ nodes which correspond to degrees of freedom: if we use the array $M$ introduced in §A-1, which is now $n_D \times n_T$ where $n_T$ is the number of elements, we have

$$v(F^j(\xi, \eta)) = \sum_{i=1}^{n_D} p_i(\xi, \eta) v_{m_{ij}}^{m_{ij}}$$

on the $j$-th element $K$, where $F^j$ denotes the mapping $\hat{K} \to K$.

In the case of isoparametric elements, the mapping $F_j$ can be evaluated explicitly from the coordinates of the nodes $a_i$ of the element:

$$F^j(\xi, \eta) = \sum_{i=1}^{n_D} a_i p_i(\xi, \eta)$$

The Jacobian matrix $\nabla F^j$ is now depending on $\xi, \eta$ (unless the element is usually "straight").

The formula (17) still holds, with a matrix $B$ now depending on $\xi, \eta$. The stiffness submatrix is now $n_D \times n_D$.

From the formula (17) it can be seen that the stiffness submatrix $A^K$ corresponding to the term $\int_K \nabla u \cdot \nabla v dx$ can be evaluated from the values of the partial derivatives of the basis functions at the integration nodes $b_\ell$ of the reference element.

Finally the data which are characteristic of the type of interpolation (or element) and necessary for the computation of the stiffness submatrices (and the corresponding terms necessary for the computation of the right hand side) are
• the coefficients $\hat{W}_{\ell}$

• the values of the basis functions at the integration nodes $\hat{b}_{\ell}$ on the reference element

• the values of their partial derivatives on the reference element

The flexibility of this system of data representation explains why, for sophisticated elements, one never uses the exact formula (22), even if the elements are straight.

### Exercise 8:

Consider the case of isoparametric Lagrange quadrilateral finite elements Q1(Fig. 9).

<figure><img src="image_16.png" /><figcaption>Fig. 9 Quadrilateral bilinear lagrange finite elements</figcaption></figure>

1. Compute the basis functions $(p_i(\xi, n))_{1 \le i \le 4}$.

2. Let $\begin{pmatrix} x_i \\ y_i \end{pmatrix}$ be the coordinate of the nodes $a_i$, $1 \le i \le 4$.

We choose the simplest integration formula: $L = 1$

$$ \hat{W}_1 = 1 \quad \hat{b}_1 = (\frac{1}{2}, \frac{1}{2}) \quad (= \text{center of } \hat{K}). $$

Compute $B = \nabla F(\hat{b}_1)$ with respect to $(x_i, y_i)$ $1 \le i \le 4$.
3. Consider the case where

$$a_1 = \begin{pmatrix} \alpha \\ \beta \end{pmatrix}, \quad a_2 = \begin{pmatrix} \alpha + h \\ \beta \end{pmatrix}, \quad a_3 = \begin{pmatrix} \alpha + h \\ \beta + h \end{pmatrix}, \quad a_4 = \begin{pmatrix} \alpha \\ \beta + h \end{pmatrix}$$

and compute the corresponding stiffness submatrix $A^K$ (associated to the Dirichlet problem as in Exercises 5.

4. Same question as in exercise 6. (All the elements are now simple squares, the side of which is h).

## C - Computation of the Triangulation

The first (and often the simplest) method to compute a triangulation of a domain is to make a special program. Of course the output of this program will be the arrays M, Z and eventually L described in §A-1.

But, for complicated domains the task can be rather tedious. Note that in any case it will be necessary to describe in a convenient way the boundary of the domain, and even that can be difficult.

We shall mention 2 methods:

### 1. Automatic Retriangulation [2] [3]

Let $T_0$ be a simple triangulation of the domain $\Omega$, it is always possible to deduce a thinner triangulation $T_1$ by the following process.

<figure><img src="image_17.png" /><figcaption>Fig. 10 Initial triangulation $T_0$</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 11 Triangulation $T_1$</figcaption></figure>

In the case of a polygonal domain, the process is very simple: you join all the middle points of the edges of the triangulation, in such a way that each triangle is divided into 4 triangles (Fig. 11).

(The Sierpinski triangle given p.105 performs this).

(b) Ω is not polygonal

<figure><img src="image_19.png" /><figcaption>Fig. 12 Retriangulation in the case of a non-polygonal domain</figcaption></figure>

In the case of a non-polygonal domain the middle point of a boundary edge is no more on the boundary Ω, so that the triangulation $T_1$ defined by the trick used in the polygonal case does not approximate accurately the boundary of Ω.
In such a case it is convenient to project each
middle point of 2 boundary nodes onto the boundary; you get
a new triangulation $T_1'$ which is more accurate (Fig. 12)
than the previous one $T_0$.

In the same way, from $T_1'$, one would deduce another
triangulation $T_2'$, and so on ...

We see that with a correct description of the
boundary of Ω which is necessary in any case for the
projection of a point onto the boundary, it is possible to
generate with simple initial data (the rough triangulation $T_0$)
a rather sophisticated triangulation.

2. A Direct Method [3]

We shall describe a method which uses as input only
a certain number of points given on the boundary.

We shall consider for simplicity only the case of
a simply connected domain, so that the initial data is a
closed polygonal line $\Gamma_0$.

By induction we assume that a polygonal line $\Gamma_n$ ($n \ge 0$)
is given, we shall construct a smaller one $\Gamma_{n+1}$, and so on ...
The algorithm will stop either when $\Gamma_n$ is reduced to a simple
triangle, or after a number of steps that you don't want to
exceed.

To construct $\Gamma_{n+1}$ you check first if there exists
2 consecutive side of $\Gamma_n$ forming on angle smaller than $\frac{\pi}{2}$.
In such a case (Fig. 13), you join the 2 opposite vertices
in order to form a triangle; you "pull out" this triangle
and get $\Gamma_{n+1}$.
<figure><img src="image_20.png" /><figcaption>Fig. 13 Case where 2 consecutive sides of $\Gamma_n$ make an angle smaller than $\pi/2$</figcaption></figure>

In case the smallest angle between 2 consecutive sides of $\Gamma_n$ is greater than $\pi/2$, create a new vertex in this angle (Fig. 14), join the new vertex to the 3 points defining the 2 consecutive vertices, and pull out the 2 triangles thus created in order to get $\Gamma_{n+1}$.

<figure><img src="image_21.png" /><figcaption>Fig. 14 Case where the smallest angle between 2 consecutive sides exceeds $\pi/2$.</figcaption></figure>
Of course what we have said about this algorithm is just
sufficient to make the reader understand the idea, but
not define it precisely. We refer the interested reader
to GEORGE [3] who originates the method.

The problem is of course that the algorithm may
fail to converge within a reasonable number of iterations.
In order to avoid this the programming has to be made
carefully.

What we can say in conclusion of this part is that
the computation of a triangulation may be the most difficult
part of the programming of a finite element problem, on a
complicated domain. This, of course, is not due to the
nature of the finite element method, rather to the complication
of the problem. (The use of a finite difference method in
such a case is very much more complicated).

D

-

# Description of 2 Methods for the Solution of Large Problems

## 1. The Frontal Method (IRONS [4])

We shall review the properties of Gaussian elimination
for the matrices with a band-structure (see Ch. 1, §3).

### Property 1:

Let $w$ denote the bandwidth of $A$, the matrix $A^{(k)}$ at the $k$-th step of the Gaussian elimination has the same bandwidth $w$. Moreover the number of collective terms to add to $A^{(k)}$ in order to get $A^{(k+1)}$ is $w^2$.

### Proof:

We assume by induction that $A^{(k)}$ has a bandwidth $w$.
We see from the formula (33) (Ch. 1) defining $A^{(k+1)}$ that
to get $a_{ij}^{(k+1)}$ we have to add to $a_{ij}^{(k)}$ a corrective term

$$ - \frac{a_{ik}^{(k)} \cdot a_{kj}^{(k)}}{a_{kk}^{(k)}}, \quad \text{for } i, j = k+1, \dots, N. $$

As a matter of fact

the non-zero corrective terms are obtained only for

$$ i, j = k+1, \dots, k+w $$

since by assumption $A^{(k)}$ has a bandwidth w.

On the other hand we see that $A^{(k+1)}$ will have the same bandwidth w as $A^{(k)}$ since none of these corrective terms are outside the band.

The Property 1 shows that at each step of the Gaussian elimination, only a small submatrix w x w of A changes.

As a matter of fact, it happens that the assembling algorithm described in section A (§2) has the same property.

**Property 2:**

Assume that the nodes of the finite element mesh are numbered in such a way that the stiffness matrix A has a bandwidth w. If we perform the assembling algorithm in N steps (where N is the number of degrees of freedom, instead of $n_T$ where $n_T$ is the number of triangle) in the following way.

At the k-th step (k = 1, ... N), assemble all the elements containing the degree of freedom k (and not already assembled in a previous step).
Then only the elements of the submatrix
(w + 1) x (w + 1)

$$
\begin{pmatrix}
a_{kk} & \cdots & a_{k,k+w} \\
& \ddots & \\
a_{k+w,k} & \cdots & a_{k+w,k+w}
\end{pmatrix}
$$

of the "current" stiffness matrix change at the $k^{th}$ step.
Moreover they are changed by the addition of a corrective
term.

The idea of the frontal method lies on the combination
of these 2 properties, and on the commutativity of the addition.
We can perform together the assembling <u>and</u> the Gaussian elimi-
nation of the stiffness matrix.

At the $k^{th}$ step of the frontal method, the (k-1) first
lines of U (and columns of L) have already been computed.
Only the submatrix

$(A_{ij})_{i,j = k, k+1, ..., k+w}$

is "active" (that is, containing non-zero elements still to
be changed).

At last, the remaining part of the matrix is still
zero, and does not yet need to be stored in the memory of the
computer.

As a matter of fact, only the "active" part needs to
be in the fast memory of the computer, that is an array of
size (w + 1) x (w + 1).
In Fig. 16-17, we show on a simple example how the frontal method works. As in Fig. 4, we denote symbolically by a small 'a' all the elements of the 4 x 4 stiffness sub-matrix of the element a; b for the element b, and so on.

In the same way, we shall denote symbolically by 1 all the corrective terms resulting from the 1st Gaussian elimination step, 2 for the 2nd etc. ... As in Fig. 4, all the terms which are in the same square representing an element of the matrix A (Fig. 17) are supposed to be added.

<figure><img src="image_22.png" /><figcaption>Fig. 16 The triangulation for Fig. 17</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 17 State of the matrix A at the 3rd step of the frontal algorithm. (The active part is the n x 6 submatrix (((a_{ij}))_{ij = 3,4,...,8}).</figcaption></figure>

For an n x n mesh, we have N = n² degrees of freedom. The bandwidth w = n, so that for large n, the number of operations for solving the problem is approximatively n⁴.

## §2. The Substructure Method

As a matter of fact the frontal method is a strategy for the choice of the order of elimination for the Gaussian elimination of a stiffness matrix. We shall consider here a completely different strategy which gives less computation for very large systems.
## 2.1 Elementary Example

Let us consider a triangulation with four square elements, and nine nodes numbered in such a way that the nodes 1 to 4 are the corner nodes (Fig. 18).

<figure><img src="image_24.png" /><figcaption>Fig. 18 Stiffness matrix in an elementary example</figcaption></figure>

We can see that, in the associated stiffness matrix A, the first four unknowns, associated to the corner nodes, are weakly related to the others. For example the Gaussian elimination of these 4 variables is easy. In fact, the Gaussian elimination of the first variable produces corrective terms (occurring in the lines and columns 5 to 9) which can be computed only from the elements of the stiffness submatrix of the element a. Moreover, these corrective terms will act only on terms containing already a contribution from the element a (in our example the lines and columns 5 to 7).
By using again the commutativity of the addition, we deduce that the matrix $A^{(5)}$, obtained from the matrix A after the elimination of the first four variables, could have been obtained directly by assembling together the 3 x 3 stiffness submatrices partially eliminated (Fig. 19).

In Fig. 19, we denote (always symbolically) by

$$ \begin{pmatrix} a' & a' & a' \\ a' & a' & a' \\ a' & a' & a' \end{pmatrix}, \begin{pmatrix} b' & b' & b' \\ b' & b' & b' \\ b' & b' & b' \end{pmatrix}, \dots $$

the 3 x 3 matrices obtained from the stiffness submatrices of the elements a, b, ..., after the Gaussian elimination of the variable relative to the corner node.

<figure><img src="image_25.png" /><figcaption>Fig. 19 Result of a direct assembling of the partially eliminated stiffness submatrices</figcaption></figure>
We shall say that ((a')) is the stiffness submatrix of the element a, after elimination of the variable number 1. Note that for the following, we can consider that ((a')) is the stiffness submatrix of the element a, if we forget the node 1 (this explains the new way of drawing the triangulation on Fig. 19).

## 2.2 Generalization

We consider now a more sophisticated example (Fig. 20).

<figure><img src="image_26.png" /><figcaption>Fig. 20 Initial triangulation</figcaption></figure>

We shall group the 16 elements (4 by 4) in four macro elements A, B, C, D (with A = a U b U c U d, ...).

We can get the 9 x 9 stiffness submatrix of the element A by assembling the four elements a, b, c, d. In the matrix A obtained, we can then eliminate the variable corresponding to the 1st node (Cf. Fig. 20), and obtain an 8 x 8 matrix A'.

We proceed in an analogous way for the macro elements B, C, D (Fig. 21).
<figure><img src="image_27.png" /><figcaption>Fig. 21 Triangulation after partial assembling and partial elimination</figcaption></figure>

We can then begin again the process by assembling the four
macro elements in a greater one E. In the corresponding
21 x 21 matrix, we eliminate again the variables corresponding
to the internal nodes (5, 6, 7, 8, 9).

We get then a 16 x 16 matrix E' corresponding to
the triangulation of Fig. 22.

<figure><img src="image_28.png" /><figcaption>Fig. 22 Triangulation after assembling
of the macro elements A'B'C'D'
and elimination of the internal
nodes</figcaption></figure>
To solve the initial system, we would then have just to solve a
linear system of simultaneous linear equations with matrix E'.

Note that there are many possible strategies for the choice
of the variables to eliminate. In the case of Fig. 21, we could
have kept the nodes 12, 5, 7, 9, 23 instead of the boundary nodes
for the last step. This strategy would have produce in fact less
computations, in this special case.

However the strategy of always eliminating the internal
nodes has the advantage to be extended immediately to a regular
mesh $2^P \times 2^P$ where p is any positive integer.

In fact, we shall perform the previous partial assembling
and elimination four elements by four elements recursively.

By induction at the $k^{th}$ step, we shall have $4^{P-k+1}$
elements, that we shall group 4 by 4, in order to get $4^{P-k}$
macro elements (Fig. 23).

<figure><img src="image_29.png" /><figcaption>Fig. 23 Transformation of the $k^{TH}$ step</figcaption></figure>

In each macro element we shall perform the elimination of the
$2^{k+1} - 1$ internal variables (Fig. 23).
Let $n = 2^p$, the number of operations of such a method is asymptotically $O(n^3)$ (compare to the $O(n^4)$ of the frontal method).

However the size of the memory necessary for the solution of the *last* system of linear equations is $4n \times 4n$ (instead of $n \times n$ for the frontal method).

The programming of such a method will then be more complicated. However all the symmetries of the structure can be more easily taken into account.

Note that the *last* linear system of simultaneous equations has only boundary variables, exactly as in the *boundary integral method* for solving an elliptic problem.

This method is called the *substructure method*. Let us mention again that it can be easily extended to higher degree elements.

For elements with internal nodes, note that it is always possible (even if you use the frontal method) to *eliminate* the internal nodes when one computes the stiffness submatrix of the element. Such a process is called "static condensation", and is different of what happens for serendipity elements.
References for Chapter 3

[1] LIONS, J.L.: These lecture notes.

[2] BOURGAT, J.F.: "Analyse Numérique de la torsion élasto-plastique", Thèse, Paris, 1970.

[3] GEORGE, J.A.: Computer implementation of the finite element method, Ph.D. Thesis, Stanford University, 1970.

[4] IRONS, B.M.: "A frontal solution program for finite element analysis", Int. J. Num. Meth. Eng 2(1970) 5-32.

[5] GEORGE, J.A.: "Block elimination on finite element systems of equations" in symposium on sparse matrices and their applications, D.J. ROSE, R.A. WILLOUGHBY, Ed., Plenum Press (1972) 101-114.
# Annex to Chapter 3

## Listing of the following subroutines

FEM1CF: FEM of degree 1 + Choleski factorization

FEM2RM: FEM of degree 2 + Relaxation method

TR1 : Construction of a triangulation for FEM2RM by refinement of a coarse triangulations.
SUBROUTINE FEM1CF(UO,F,Q,ME,IRA,NG,A,C2,NT,NS,MUD,MM,KSW,U,BF,CQ)

COMPUTES THE SOLUTION OF -DELTA U= F, IN OMEGA, U ON THE BOUNDARY = UO
BY FINITE ELEMENTS OF ORDER 1 WITH CHOLESKI FACTORISATION OF THE MATRIX.

INPUT
    MM      NUMBER OF INTERIOR NODES
    NS      NUMBER OF NODES
    NT      NUMBER OF TRIANGLES
    Q(2,NS) COORDINATES OF THE NODES
    ME(3,NT) ARRAY OF THE NUMBERS OF THE VERTICES OF THE TRIANGLES
    IRA(NS) NUMBERS OF THE NODES WHEN THE BOUNDARY NODES ARE NOT COUNTED
    NG(NS) =1 IF THE NODE IS ON THE BOUNDARY =0 OTHERWISE
    MUD     BANDWIDTH OF THE MATRIX
    KSW=1    F IS GIVEN ON THE NODES
             =4 F IS PIECE-WISE CONSTANT ON THE TRIANGLES
             =6 SUM ON OMEGA OF F*BASIS FUNCTION IS GIVEN AS INPUT
    F(NS OR NT) ACCORDING TO KSW
    IF KSW=2,3,5 C2,CQ,A FACTORIZED ARE INPUT AND F IS AS ABOVE
    UO(NS) VALUES OF U ON THE BOUNDARY

IF NG=0 UO AT THIS POINT IS NOT IMPORTANT (ONLY USED AS WORKING SPACE)
OUTPUT
    U(NS) SOLUTION
    BF(NS) SUM ON OMEGA OF F*BASIS FUNCTION
    C2(NT) 2*AREA OF TRIANGLES
    CQ(6,NT) DERIVATIVES OF BASIS FUNCTIONS ON EACH TRIANGLE
    A((MUD+1)*(2*M-MUD)/2) MATRIX FACTORISED

WORKING ARRAYS DP(6),BF(NS),UO(NS)
ON OUTPUT F AND UO ARE RESTITUTED.

COMMENTS IN FRENCH

DIMENSION UO(1),F(1) ,ME(3,1),IRA(1),NG(1),A(1),CQ(6,1) ,C2(1)
DIMENSION U(1),BF(1),DP(6),Q(2,1)
X24=1./24.
DP(1)=-1.D0
DP(2)=-1.D0
DP(3)=1.D0
DP(4)=0.D0
DP(5)=0.D0
DP(6)=1.D0
X6=1./6.

FACTORISATION DE LA MATRICD

IF(KSW.EQ.3.OR.KSW.EQ.2.OR.KSW.EQ.5)GO TO 2
CALL ASTOF(Q,ME,IRA,A,BF,DP,CQ,C2,NG,NT,MUD,MM)
CALL DMCHB(I),A,MM,1,MUD,0,1,E-15,IER)
PRINT 100,IER
100 FORMAT(I3)
CALCUL DU SECOND MEMBRE
2 CONTINUE
DO 1 I=1,NS
BF(I)=0.D0
1 CONTINUE
DO 10 JK=1,NT
KK=-2
DO 10 KC=1,3
KK=KK+2
JJ=IRA(ME(KC,JK))
IF(JJ.EQ.0)GO TO 22
K=-2
DO 13 KB=1,3
K=K+2
J=ME(KB,JK)
IF(NG(J).NE.1)GO TO 12
UNG=UO(J)
IF(UNG.EQ.0,DO)GO TO 12
DO 11 L=1,2

11 BF(JJ)=BF(JJ)*UNG*CQ(K+L,JK)*CQ(KK+L,JK)*C2(JK)/2.
12 CONTINUE
13 CONTINUE
22 CONTINUE
10 CONTINUE
15 CONTINUE
DO 14 I=1,NS
14 U(I)=0.D0
IF(KSW.GE.5)GO TO 33

IF(KSW.GE.3)GO TO 30

DO 3 J=1,NT

C2J=C2(J)

ME1=ME(1.J)

ME2=ME(2.J)

ME3=ME(3.J)

U(ME1)=U(ME1)•C2J*x24*(2.*F(ME1)+F(ME2)+F(ME3))

U(ME2)=U(ME2)•C2J*x24*(F(ME1)+2.*F(ME2)+F(ME3))

U(ME3)=U(ME3)•C2J*x24*(F(ME1)+F(ME2)+2.*F(ME3))

3 CONTINUE

GO TO 31

30 DO 32 J=1,NT

DO 32 K=1,3

KJ=ME(K,J)

32 U(KJ)=U(KJ)•F(J)*C2(J)*X6

GO TO 31

33 DO 34 I=1,NS

34 U(I)=F(I)

31 CONTINUE

DO 4 I=1,NS

IF(IRA(I).EQ.0)GO TO 8

S=U(I)

BF(IRA(I))=S•BF(IRA(I))

8 CONTINUE

4 CONTINUE

CALL DMCHB(BF.A,MM,1,MUD,-1,1,E-15,IER)

IF(IER.NE.0)PRINT 100,IER

DO 9 I=1,NS

9 IF(IRA(I).NE.0)UO(I)=BF(IRA(I))

DO 24 I=1,NS

BF(I)=U(I)

U(I)=UO(I)

IF(IRA(I).NE.0)UO(I)=0.00

24 CONTINUE

RETURN

END
SUBROUTINE DPQH

PURPOSE: TO COMPUTE THE STIFFNESS SUBMATRIX B4 OF ONE ELEMENT K, RELATIVE TO THE BILINEAR FORM:

DESCRIPTION OF PARAMETERS

INPUT: DP(2,3) VALUES OF THE DERIVATIVES OF THE BASIS FUNCTIONS ON THE REFERENCE ELEMENT
        X(2,3) COORDINATES OF THE VERTICES OF THE ELEMENT K

OUTPUT: DQ(2,3) VALUES OF THE DERIVATIVES OF THE BASIS FUNCTIONS ON THE ELEMENT
         D2 2 TIMES THE AREA OF THE ELEMENT
         B4(6) THE STIFFNESS SUBMATRIX OF THE ELEMENT UNDER SYMMETRIC COMPRESSED FORM (CONTAINS B11,B21,B22,B31,B32,B33)
                 ( TO BE MULTIPLIED BY THE AREA OF THE ELEMENT)
         DF1 LENGTH OF THE SIDE 1-2.

SUBROUTINE DPQH(DP,X,DQ,D2,R4,DF1)

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

DIMENSION DP(1),X(1),DQ(1),B4(1),G1(4),G3(2,2),SJB(4)

DATA SJB/1.D0,-1.D0,-1.D0,1.D0/

FORMATION DU JACOBIEN

IL=0
D0 2 I=1,2
D0 2 L=1,2
S=0.D0
KC=-2
DQ 3 K=1,3
KC=KC+2

3 S=S+DP(I+KC)*X(L+KC)

IL=IL+1

2 G1(IL)=S

FORMATION DE G3 MAT. DES COFACTEURS DE G1

IL=0
D0 4 I=1,2
D2=0.D0
D0 4 L=1,2
IL=IL+1
G3(L,I)=G1(5-IL)*SJB(IL)
D2=D2+G3(L,I)*G1(IL)
CALCUL DE DQ ET DE B3
KC=-2
D0 7 K=1.3
KC=KC+2
D0 5 I=1,2
S=0.D0
D0 6 L=1,2

S=S+G3(I,L)*DP(L+KC)

DQ(I+KC)=S/D2

CONTINUE

CALCUL DE LA SOUS MATRICE B4

IK12=0
<1C=-2
D0 18 K1=1,3
K1C=K1C+2
<2C=-2
D0 18 K2=1,3
K2C=K2C+2
IF(K2-K1)11,11,18
IK12=IK12+1
S=0.D0
D0 9 IK=1,2

3 S=S+DQ(IK+K1C)*DQ(IK+K2C)

B4(IK12)=S

18 CONTINUE

AIRE DE LA FACE 1

DF1=0.D0

D0 1 L=1,2

DF1=DF1*G3(L,1)*G3(L,1)

DF1=DSQRT(DF1)

RETURN

END
SUBROUTINE ASTOF

PURPOSE: ASSEMBLE THE STIFFNESS MATRIX A UNDER COMPRESSED SYMMETRIC FORM, READY TO BE USED BY IBM SSP SUBROUTINE DMCHB.

PARAMETERS: SEE FEM1CF.

USES DPQH.

SUBROUTINE ASTOF(Q,ME,IRA,A,BF,DP,CQ,C2,NG,NT,MUD,M)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION MF(3),G4(12),GX(6),Q(1),ME(1),DP(1),A(1),IRA(1),BF
1(1),CQ(1),C2(1),NG(1),NGL(3),UG(1)
MA=(MUD*1)@(2*M-MUD)/2
PRINT 805,M,MUD

805 FORMAT(' M=',I4,' MUD=',I4)
PRINT 874,MA

874 FORMAT(' MA= ',I5)
D0 40 I=1,M

40 BF(I)=0.00
D0 2 LA=1,MA

2 A(LA)=0.00

7 NUD=0
JQ=-5
D0 6 J=1,NT
JQ=JQ+6
JH=3*(J-1)
KI=0

D0 11 K=1,3
MEK=ME(K+JB)
NGL(K)=NG(MEK)
MF(K)=IRA(MEK)
D0 11 IB=1,2
KI=KI+1
MFL=IB+2*(MEK-1)

11 GX(KI)=Q(MFL)
CALL DPQH(DP,GX,CQ(JQ),C2(J),G4,GF1)
KI=0
D0 12 KI=1,3
D0 20 L=1,KI
MJ=MF(KI)
KL=KL+1
MI=MF(L)
IF(MI.LE.0.OR.MJ.LE.0)GO TO 20
IF(MJ-MI)27,28,28

27 MS=MJ
MJ=MI
MI=MS

28 IF((MJ-MI).GT.MUD)PRINT 872,J
NUD=MAX0(NUD,MJ-MI)
MS=MAX0(MI-1-M*MUD,0)
MS=MS*(MS+1)/2
LA=(MI-1)*MUD+MJ-MS
IF((LA#(1+MA-LA)).LE.0)PRINT 872,J
A(LA)=A(LA)+G4(KL)*C2(J)/2.D0

20 CONTINUE

12 CONTINUE

872 FORMAT(' ERREUR D INDICE J= ',I4)

6 CONTINUE
IF(MUD.EQ.NUD)RETURN
NA=(NUD+1)*(2*M-MUD)/2
PRINT 873,NUD,NA

873 FORMAT(' ON PEUT PRENDRE MUD= ',I4,' ET MA= ',I8)
RETURN
END
C

SUBROUTINE FEM2RM(NSEGVM,IOPT,NDIV,EPSOM,EPSREL,NT,NS,NSI,C,NU,AXX
1,AXY,AXX,AYX,AO,NAP,F,ISW,U,NSEGV,NUSEGV,ISEG,DIA,A,B,AIRE,NSEG,
1 NSEGI,JSEG)

C

COMPUTES THE SOLUTION OF

C

-DX(AXXDX(U)+AYXDY(U))-DY(AXYDX(U)+AYYDY(U))=0 IN OMEGA

C

U ON BOUNDARY =0.

C

BY LAGRANGIAN TRIANGULAR FINITE ELEMENTS OF DEGRE 2

C

THE LINEAR SYSTEM IS SOLVED BY THE RELAXATION METHOD

C

INPUT

C

EPSOM PRECISION PARAMETER . A GOOD VALUE IS 0.1 (FOR COMPUTATION OF OMEGA

C

EPSREL PRECISION PARAMETER TAKE= 10**-5 (STOP RULE FOR RELAXATION ALGO)

C

AXX(NT),AXY(NT),AYX(NT),AYY(NT),AO(NT) (GIVEN FOR REFINED TRIANG.)

C

COEFFICIENTS OF PDE. ASSUMED PIECEWISE CONSTANTS ON TRIANGLES

C

NAP(NS-NSI) =0 IF NODES IS BOUNDARY-DIRICHLET ,=1 OTHERWISE.

C

F(NS) RIGHT MEMBER OF PDE. PIECEWISE LINEAR.

C

ISW=2 PRINTING OF SOLUTION =1 NO PRINTINGS.

C

GIVEN AS OUTPUT OF TRI NSEGVM,NT,NS,NSI,IOPT,C,NU,NDIV

C

NUSEGV(NS),NUSEGV(8,NS),ISEG(2,NS),JSEG(2,NT),DIA(NS),A(2,NT),B(NS),AIR(NT)

C

OUTPUT

C

U(NS) SOLUTION

C

DIMENSION C(2,1),NU(3,1),AXX(1),AXY(1),AYX(1),AYY(1),AO(1),F(1)

C

DIMENSION U(1),NSEGV(1),NUSEGV(8,1),ISEG(2,1),DIA(1),A(2,1),B(1),

C

1 AIRE(1)

C

DIMENSION NAP(1) ,JSEG(2,1)

C

MATRICE

C

********

C

CALL AIRET(NT,C,NU,AIRE)

C

CALL MAT1(NT,NS,NSEG,NSEGI,C,NU,ISEG,JSEG,AIRE,AXX,AYY,AXY,AYX,
1 AO,DIA,A)

C

CONDITIONS AUX LIMITES

C

*********************

C

SECOND MEMBRE

C

*********************

C

CALL FW1(NT,NS,NU,AIRE,F,1,B,1)

C

RESOLUTION PAR RELAXATION

C

*********************

C

CALL OMEGA1(NS,NSI,NSEGVM,NSEGV,NUSEGV,ISEG,NAP,DIA,A,U,OMEGA,
1 EPSOM,NITOM)

C

DO 40 I=1,NS

C

U(I)=0.

40

CONTINUE

C

CALL RELAX1(NS,NSI,NSEGVM,NSEGV,NUSEGV,ISEG,NAP,DIA,A,R,U,OMEGA,
1 EPSREL,NIT)

C

IMPRESSION DES RESULTATS

C

*********************

C

IF(ISW,NE,2)RETURN

C

PRINT 70,NIT

C

FORMAT(1H1,1X,'SOLUTION APPROCHEE',5X,'ITERATIONS DE RELAXATION',
114,7)

C

PRINT 61,(U(I),I=1,NS)

61

FORMAT(20F6.3)

RETURN

C

END

C

SUBROUTINE AIRET(NT,C,NU,AIRE)

C

*********************

C

AIRE DES TRIANGLES D UNE TRIANGULATION

C

DIMENSION C(2,1),NU(3,1),AIRE(1)

C

DO 10 J=1,NT

C

I1=NU(1,J)

C

I2=NU(2,J)

C

I3=NU(3,J)

C

C11=C(1,I1)

C

C12=C(1,I2)

C

C13=C(1,I3)

C

C21=C(2,I1)

C

C22=C(2,I2)

C

C23=C(2,I3)

C

AIRE(J)=0.5*(C11*C22-C12*C21+C12*C23-C13*C22+C13*C21-C11*C23)

C

CONTINUE

RETURN

C

END

10
SUBROUTINE MAT1(NT,NS,NSEG,NSEGI,C,NU,ISEG,JSEG,AIRE,AXX,AYY,AXY,
1AYX,A0,DIA,A)

TERMES NON NULS DE LA MATRICE D UN PB ELLIPTIQUE D ORDRE 2

DISCRETISE PAR ELEMENTS FINIS D ORDRE 1 SUR UNE TRIANGULATION

AUTOMATIQUE DE TYPE 1. STOCKAGE PAR DIAGONALE ET SEGMENTS.

DIMENSION C(2,1),NU(3,1),ISFG(2,1),JSEG(2,1),AIRE(1),AXX(1),
1AYY(1),AXY(1),AYX(1),A0(1),DIA(1),A(2,1)

DIMENSION IP1(3),KP1(2)

IP1(1)=2

IP1(2)=3

IP1(3)=1

KP1(1)=2

KP1(2)=1

DO 10 I=1,NS

DIA(I)=0.

CONTINUE

DO 20 J=1,NT

DO 20 K=1,3

I=NU(K,J)

K1=IP1(K)

K2=IP1(K1)

I1=NU(K1,J)

I2=NU(K2,J)

AA=C(2,I1)-C(2,I2)

BR=C(1,I2)-C(1,I1)

DIA(I)=DIA(I)+0.25*(AXX(J)*AA*AA+AYY(J)*BR*BR+(AXY(J)+AYX(J))*AA*
1HR)/AIRE(J)+AIRE(J)*A0(J)/6.

CONTINUE

DO 30 L=1,NSEG

A(1,L)=0.

A(2,L)=0.

CONTINUE

DO 40 L=1,NSEG

DO 50 K=1,2

IF(L.GT.NSEGI.AND.K.EQ.2) GOTO 50

I=ISEG(K,L)

I1=ISEG(KP1(K),L)

J=JSEG(K,L)

I11=NU(1,J)

IF(III.NE.I) GOTO 51

I2=NU(3,J)

GOTO 53

IF(III.NE.I1) GOTO 52

I2=NU(2,J)

GOTO 53

I2=III

AA1=C(2,I1)-C(2,I2)

AA2=C(2,I2)-C(2,I1)

BR1=C(1,I2)-C(1,I1)

BR2=C(1,I1)-C(1,I2)

CXY=AA2*BR1

CYX=AA1*BR2

COE=0.25/AIRE(J)

TT=COE*(AXX(J)*AA1*AA2+AYY(J)*BR1*BR2)+AIRE(J)*A0(J)/12.

TXY=COE*(CXY*AXY(J)+CXY*AYX(J))

TYX=COE*(CYX*AXY(J)+CXY*AYX(J))

IF(K.EQ.2) GOTO 60

A(1,L)=TT+TXY

A(2,L)=TT+TYX

GOTO 50

A(1,L)=A(1,L)+TT+TYX

A(2,L)=A(2,L)+TT+TXY

CONTINUE

CONTINUE

RETURN

END

C
C
C
C
C

C

10

20

30

51

52

53

60

50

40
C SUBROUTINE OMEGA1(NS,NSI,NSEGVM,NSEGV,NUSEGV,ISEG,NAP,DIA,A,X,
1OMEGA,EPSOM,NIT)

C RECHERCHE DU PARAMETRE OPTIMAL DE RELAXATION POUR UNE DISCRETISATION
C PAR ELEMENTS FINIS D ORDRE 1

C DIMENSION NSEGV(1),NUSEGV(NSEGVM,1),ISEG(2,1),NAP(1),DIA(1),A(2,1)
1,X(1)

C

DO 10 I=1,NS

IF(I.GT,NSI) GOTO 11

12

X(I)=1.

GOTO 10

11

IF(NAP(I-NSI).EQ.1) GOTO 12

X(I)=0.

10

CONTINUE

NIT=0

15

NIT=NIT+1

ROMA=-1000.

ROMI=1000

DO 20 I=1,NS

21

IF(I,LE,NSI) GOTO 21

IF(NAP(I-NSI).EQ.0) GOTO 20

XI=X(I)

S=0.

NSEGVI=NSEGV(I)

DO 30 K=1,NSEGVI

L=NUSEGV(K,I)

I1=ISEG(1,L)

I2=ISEG(2,L)

IF(I.EQ.12) GOTO 31

S=S+A(1,L)*X(12)

GOTO 30

31

S=S+A(2,L)*X(11)

30

CONTINUE

X(I)=-S/DIA(I)

QUO=X(I)/XI

ROMA=AMAX1(ROMA,QUO)

ROMI=AMIN1(ROMI,QUO)

20

CONTINUE

ROMA=AMIN1(ROMA,0.99999)

OMI=2./ (1. *SQRT(1.-ROMI))

OMA=2./ (1. *SQRT(1.-ROMA))

IF((OMA-OMI).GT.EPSOM.AND.NIT.LT.100) GOTO 15

OMEGA=(OMI+OMA)/2.

PRINT 40,OMEGA,NIT,EPSOM

40

FORMAT(1H1./1X,1'PARAMETRE OPTIMAL DE RELAXATION OMEGA =',F8.4,/,
11X,'NOMBRE D'ITERATIONS POUR OBTENIR OMEGA',15.5X,'AVEC LE TEST D

2ARRET EPSOM =',F8.4,/)

RETURN

END

C

SUBROUTINE PTRI(X,Y,NT,C,NU,NUTRI)

C

*******************

C

DETERMINATION D UN TRIANGLE CONTENANT UN POINT DONNE

C

DIMENSION C(2,1),NU(3,1)

C

NUTRI=-1

EPS=0.

9

DO 10 J=1,NT

I1=NU(1,J)

I2=NU(2,J)

I3=NU(3,J)

X1=C(1,I1)-X

X2=C(1,I2)-X

X3=C(1,I3)-X

Y1=C(2,I1)-Y

Y2=C(2,I2)-Y

Y3=C(2,I3)-Y

PV1=X1*Y2-X2*Y1

IF(PV1.LT.EPS) GOTO 10

PV2=X2*Y3-X3*Y2

IF(PV2.LT.EPS) GOTO 10

PV3=X3*Y1-X1*Y3

IF(PV3.LT.EPS) GOTO 10

NUTRI=J

GOTO 20

10

CONTINUE

IF(NUTRI.EQ.0) GOTO 20

NUTRI=0

EPS=-0.0001*SQRT((C(1,I1)-C(1,12))*2+(C(2,I1)-C(2,12))*2)

GOTO 9

20

RETURN

END
SUBROUTINE RELAX1(NS,NSI,NSFGVM,NSEGV,NUSEGV,ISEG,NAP,DIA,A,B,X,
1OMEGA,EPSREL,NIT)

C
C
C
C

RESOLUTION PAR LA METHODE DE RELAXATION D UN SYSTEME LINEAIRE
ORTENU PAR ELEMENTS FINIS D ORDRE 1

DIMENSION NSEGV(1),NUSEGV(NSEGVM,1),ISEG(2,1),NAP(1),DIA(1),
1A(2,1),B(1),X(1)

C

NIT=0

10

NIT=NIT+1

S1=0.

S2=0.

DO 20 I=1,NS

IF (I,LE,NSI) GOTO 21

IF (NAP(I-NSI),EQ,0) GOTO 20

21

XI=X(I)

S=0.

NSEGV1=NSEGV(I)

DO 30 K=1,NSEGVI

L=NUSEGV(K,I)

I1=ISEG(1,L)

I2=ISEG(2,L)

IF (I,EQ,12) GOTO 31

S=S+A(1,L)*X(I2)

GOTO 30

31

S=S+A(2,L)*X(I1)

30

CONTINUE

X(I)=(-S+H(I))/(DIA(I))

X(I)=XI+OMEGA*(X(I)-XI)

S1=S1+ABS(XI-X(I))

S2=S2+ABS(X(I))

20

CONTINUE

DRM=S1/S2

IF (DRM,GT,EPSREL) GOTO 10

RETURN

END
SUBROUTINE TRI(NDIV,NT,NS,NSI,NSEG,NSEGI,C,NU,NSEGVM,NSEGV, 1,NUSEGV,ISEG,JSEG,IOPT)

AN INITIAL COARSE TRIANGULATION IS GIVEN AND THE PROGRAM SUBDIVID IT.
TRIANGULATION

WILL BE DIVIDED INTO 4**NDIV TRIANGLES

INPUT

NDIV NUMBER OF REFINEMENT OF THE INITIAL TRIANGULATION. EACH TRIANGLE

IOPT =1 MINIMUM OUTPUT FOR THE FINAL TRIANGULATION

2 PRINT THE COORDINATES OF THE NODES

3. PRINT COORDINATES AND NUMBERS OF VERTICES OF TRIANGLES

NSEGVM MAXIMUM NUMBER OF SIDES HAVING A COMMON VERTEX. GREATER THAN 6 PLEASE

NT NUMBER OF TRIANGLES OF INITIAL TRIANGULATION

NS NUMBER OF NODES OF INITIAL TRIANGULATION

NSI NUMBER OF INTERIOR NODES OF INITIAL TRIANGULATION

C(2,NS) COORDINATES OF NODES OF INITIAL TRIANGULATION

BOUNDARY NODES MUST BE NUMBERED LAST

NU(3,NT) ARRAY OF NUMBERS OF VERTICES OF TRIANGLESS OF INITIAL TRI.

OUTPUT /RELATIVE TO THE REFINED TRIANGULATION)

NT,NS,NSI,NSEG,NSEGI,C,NU,NSEGV,NUSEGV,ISEG,JSEG

SUBROUTINE GERANT SEGV,DIVT,RESTRI.POUR UNE TRIANGULATION

DIMENSION C(2,1),NU(3,1),NSEGV(1),NUSEGV(NSEGVM,1),ISEG(2,1),
1,JSEG(2,1)

DO 10 K=1,NDIV

CALL SEGV(NT,NS,NSI,NSEG,NSEGI,NU,NSEGVM,NSEGV,NUSEGV,ISEG,JSEG)

CALL DIVT(NT,NS,NSI,C,NU,ISEG,JSEG)

CONTINUE

CALL SEGV(NT,NS,NSI,NSEG,NSEGI,NU,NSEGVM,NSEGV,NUSEGV,ISEG,JSEG)

CALL RESTRI(NT,NS,NSI,NSEG,NSEGI,C,NU,IOPT)

RETURN

END

SUBROUTINE RESTRI(NT,NS,NSI,NSEG,NSEGI,C,NU,IOPT)

**************************************************

RESULTATS DE LA TRIANGULATION

DIMENSION C(2,1),NU(3,1)

PRINT 10

FORMAT (1H1,'TRIANGULATION FINALE',/,.1X,20('*'),/)

PRINT 20,NT,NS,NSI,NSEG,NSEGI

FORMAT (1X,'NOMBRE DE TRIANGLES',16,/,.1X,'NOMBRE DE NOEUDS',16,/,
11X,'NOMBRE DE NOEUDS INTERIEURS',16,/,.1X,'NOMBRE DE SEGMENTS',16,
2/,.1X,'NOMBRE DE SEGMENTS INTERIEURS',16/)

IF (IOPT.EQ.1) GOTO 70

PRINT 30

FORMAT (1H1,'COORDONNEES DES NOEUDS',/,.1X,22('*'),/)

PRINT 40,(I,C(1,I),C(2,I),I=1,NS)

FORMAT (5(15,F8.4,F8.4,4X))

IF (IOPT.EQ.2) GOTO 70

PRINT 50

FORMAT (1H1,'NUMEROS DES SOMMETS DES TRIANGLES',/,.1X,33('*'),/)

PRINT 60,(J,(NU(K,J),K=1,3),J=1,NT)

FORMAT (5(15,2X,315,3X))

RETURN

END

SUBROUTINE DIVT(NT,NS,NSI,C,NU,ISEG,JSEG)

**************************************************

DIVISION D UNE TRIANGULATION

DIMENSION C(2,1),NU(3,1),ISEG(2,1),JSEG(2,1)

DIMENSION IP1(3)

IP1(1)=2

IP1(2)=3

IP1(3)=1

NO=(3*(NT+NS)-NSI)/2

NI0=(3*(NT+NSI)-NS)/2

NF=NS-NSI

NSEG=(3*NT+NF)/2

NSEGI=NSEG-NF

ID=NI0-NSI

DO 5 I=1,NF

II=NF+1-I

C(1,NI0+II)=C(1,NSI+II)

C(2,NI0+II)=C(2,NSI+II)

CONTINUE

ON STOCKE L'ANCIEN NUI(1.J) DANS ISEG(1.J>NSEG)

5
10 CONTINUE
IF(NSEGI.EQ.0) GOTO 101
DO 100 L=1,NSEGI
I1=ISEG(1.L)
I2=ISEG(2.L)
J1=JSFG(1.L)
J2=JSFG(2.L)
I11=ISEG(1,J1+NSEG)
IF(III.NE.I1) GOTO 11
K1=1

11 GOTO 20
IF(III.NE.I2) GOTO 12
K1=3

12 GOTO 20
K1=2
20 I11=ISEG(1,J2+NSEG)
IF(III.NE.I1) GOTO 21
K2=3

21 GOTO 30
IF(III.NE.I2) GOTO 22
K2=1

22 GOTO 30
K2=2

30 KK1=IP1(K1)
KK2=IP1(K2)
J11=4*(J1-1)+K1
J12=4*(J1-1)+KK1
J13=4*J1
J21=4*(J2-1)+K2
J22=4*(J2-1)+KK2
J23=4*J2
IF(I1.GT.NSI) I1=I1+ID
NU(K1,J11)=I1
NU(KK2,J22)=I1
IF(I2.GT.NSI) I2=I2+ID
NU(K2,J21)=I2
NU(KK1,J12)=I2
I3=NSI+L
NU(KK1,J11)=I3
NU(K1,J12)=I3
NU(K1,J13)=I3
NU(KK2,J21)=I3
NU(K2,J22)=I3
NU(K2,J23)=I3
C(1,I3)=0.5*(C(1,11)+C(1,12))
C(2,I3)=0.5*(C(2,11)+C(2,12))

100 CONTINUE
NSGIPI=NSEGI+1
101 DO 200 L=NSGIPI,NSEG
I1=ISEG(1.L)
I2=ISFG(2.L)
J=JSFG(1.L)
I11=ISEG(1,J+NSEG)
IF(III.NE.I1) GOTO 110
K=1
GOTO 120
IF(III.NE.I2) GOTO 111
K=3

110 GOTO 120
K=2
120 KK=IP1(K)
J1=4*(J-1)+K
J2=4*(J-1)+KK
J3=4*J
I1=I1+ID
I2=I2+ID
NU(K,J1)=I1
NU(KK,J2)=I2
I3=NI0+NF+L-NSEGI
NU(KK,J1)=I3
NU(K,J2)=I3
NU(K,J3)=I3

200 C(1,I3)=0.5*(C(1,11)+C(1,12))
C(2,I3)=0.5*(C(2,11)+C(2,12))
CONTINUE
NT=4*NT
NS=NO
NSI=NI0
RETURN
END

200 CONTINUE
NT=4*NT
NS=NO
NSI=NI0
RETURN
END

REMARK : NN being the number of sides of the triangulation, you must declare iSEG(2,NN), jSEG(2,NN) and NSEGv(NS) and NUSEGV(NSEGVM,NS) (see example page 108)

- The subroutine SEGV is required (see P.107)
C

SUBROUTINE SEGV(NT,NS,NSI,NSEG,NSEGI,NU,NSEGVM,NSEGV,NUSEGV, ISEG,JSEG)

C

SEGMENTS VOISINS D UN SOMMET

C

DIMENSION NU(3,1),NSEGV(1),NUSEGV(NSEGVM,1),ISEG(2,1),JSEG(2,1)
DIMENSION IP1(3)

C

IP1(1)=2
IP1(2)=3
IP1(3)=1
DO 10 I=1,NS
NSEGV(I)=0
DO 10 K=1,NSEGVM
NUSEGV(K,I)=0
CONTINUE
NSEG=NT+NS-1
NSEGI=NSEG-(NS-NST)
NUSI=0
NUSF=NSEGI
DO 20 L=1,NSEG
DO 20 K=1,2
ISEG(K,L)=0
JSEG(K,L)=0

20

CONTINUE
DO 30 J=1,NT
DO 30 K=1,3
K1=IP1(K)
I1=NU(K,J)
I2=NU(K1,J)
IF(I1.GT.NSI.AND.I2.GT.NSI) GOTO 31
IF(I1.GT.I2) GOTO 30
NUSI=NUSI+1
NUS=NUSI
GOTO 32

27

DO 29 JJ=1,NT
I11=NU(1,JJ)
I12=NU(2,JJ)
I13=NU(3,JJ)
IF(III.NE.I1.AND.II2.NE.I1.AND.II3.NE.I1) GOTO 29
IF(III.NE.I2.AND.II2.NF.I2.AND.II3.NE.I2) GOTO 29
IF(JJ.EQ.J) GOTO 29
GOTO 27

29

CONTINUE

33

NUSF=NUSF+1

NUS=NUSF

32

NSEGV(I1)=NSEGV(I1)+1

NSEGV(NSEGV(I1),I1)=NUS

ISEG(1,NUS)=I1

ISEG(2,NUS)=I2

NSEGV(I2)=NSEGV(I2)+1

NUSEGV(NSEGV(I2),I2)=NUS

30

CONTINUE

DO 40 J=1,NT

DO 40 K=1,3

K1=IP1(K)
I=NU(K,J)

I1=NU(K1,J)

NSEGVI=NSEGV(I)

DO 40 KK=1,NSEGVI

L=NUSEGV(KK,I)

IF(ISEG(2,L).EQ.I1) GOTO 50

IF(ISEG(1.L).EQ.I1) GOTO 60

GOTO 40

50

JSEG(1,L)=J

GOTO 40

60

JSEG(2,L)=J

40

CONTINUE

NSGIPI=NSEGI+1

DO 70 L=NSGIPI,NSEG

JSEG(2,L)=JSEG(1,L)

70

CONTINUE

RETURN

END
## DONNEES POUR LA TRIANGULATION INITIALE

DIMENSION C(2,1),NU(3,1)

READ 10,NT,NS,NSI

FORMAT (315)

PRINT 11,NT,NS,NSI

FORMAT (1H1,'TRIANGULATION INITIALE',/.1X,22(1*'),//.1X,'NT=',15,

15X,'NS=',15.5X,'NSI=',15./)

READ 20,((C(K,I),K=1,2),I=1,NS)

FORMAT (8F8.4)

PRINT 21,(I,C(1,I),C(2,I),I=1,NS)

FORMAT (/,1X,'COORDONNEES DES NOEUDS',/,(15.6X,F9.4,2X,F9.4))

READ 30,((NU(K,J),K=1,3),J=1,NT)

FORMAT (1B14)

PRINT 31,(J,(NU(K,J),K=1,3),J=1,NT)

FORMAT (/,1X,'NUMEROS DES SOMMETS',/,(15.6X,14.2X,14.2X,14))

RETURN

END

SUBROUTINE VALF1C(NT,C,NU,XP,YP,NUTRI,COE1,COE2,COE3)

*********************

CALCULS DE COEFFICIENTS POUR L'UTILISATION DE VALF1R

DIMENSION C(2,1),NU(3,1)

CALL PTRI(XP,YP,NT,C,NU,NUTRI)

IF(NUTRI,NE,0) GOTO 10

COE1=0.

COE2=0.

COE3=0.

GOTO 20

I1=NU(1,NUTRI)

I2=NU(2,NUTRI)

I3=NU(3,NUTRI)

X1=C(1,I1)

X2=C(1,I2)

X3=C(1,I3)

Y1=C(2,I1)

Y2=C(2,I2)

Y3=C(2,I3)

A1=Y2-Y3

A2=Y3-Y1

A3=Y1-Y2

B1=X3-X2

B2=X1-X3

B3=X2-X1

C1=X2*Y3-X3*Y2

C2=X3*Y1-X1*Y3

C3=X1*Y2-X2*Y1

AIRE=0.5*(C1+C2+C3)

TT=0.5/AIRE

COE1=TT*(A1*XP+B1*YP+C1)

COE2=TT*(A2*XP+B2*YP+C2)

COE3=TT*(A3*XP+B3*YP+C3)

RETURN

END

SUBROUTINE VALF1R(NU,F,NVAL,NUTRI,COE,VAL)

*********************

VALEURS D'UNE FONCTION CONTINUE , LINEAIRE PAR MORCEAUX SUR UNE TRIANGULATION , EN DES POINTS DONNES . CETTE SUBROUTINE EST CONCUE POUR LES APPELS REPETITIFS .

DIMENSION NU(3,1),F(1),NUTRI(1),COE(3,1),VAL(1)

DO 10 K=1,NVAL

J=NUTRI(K)

IF(J,NE,0) GOTO 5

VAL(K)=0.

GOTO 20

I1=NU(1,J)

I2=NU(2,J)

I3=NU(3,J)

VAL(K)=COE(1,K)*F(I1)+COE(2,K)*F(I2)+COE(3,K)*F(I3)

CONTINUE

RETURN

END

/*

//60.SYSIN DD #

8 9 1

0.5 0.5 0.0 1.0 0.0 0.5 0.0 0.0

0.5 0.0 1.0 0.0 1.0 0.5 1.0 1.0

0.5 1.0 3 4 1 4 5 1 5 6 1 6 7 1 7 8

1 2 3 1 3 4 1 4 5 1 5 6 1 6 7 1 7 8

1 8 9 1 9 2

/*

/\*
Chapter 4

Examples of Applications

1. A linear parabolic problem

The temperature $\phi$ on a plate $\Omega$ is governed by the
parabolic equation

$$
(4.1) \quad f = \frac{\partial \phi}{\partial t} - \Delta \phi \text{ in } \Omega \times ]0, T[; \phi|_{\Gamma_1} = \phi, \frac{\partial \phi}{\partial n}|_{\Gamma_2} = g, t \in ]0, T[
$$

$$
(4.2) \quad \phi(0, x) = 0 \quad \forall x \in \Omega
$$

where $f, \phi_0, g$ are functions which represent the colorific
exchanges with the surrounding of the plate; $t$ is the time.

Let ] 0, T [ be divided into N equal interval of size $\Delta t$

Let $\phi^n = \phi(n \Delta t)$, $n = 0, 1, \dots, N$

In order to simulate (4.1) numerically we shall discretize $\frac{\partial \phi}{\partial t}$ by
by finite differences and $\Delta \phi$ by finite elements. A fully explicit
scheme would be

$$
(4.2) \quad \frac{\phi^{n+1} - \phi^n}{\Delta t} - \Delta\phi^n = f
$$

In weak form, with the notation of Chapter II it becomes

$$
(4.3) \quad \int_{\Omega} \phi^{n+1} w_i dx = \int_{\Omega} \phi^n w_i dx - \Delta t \langle \phi^n, w_i \rangle + \Delta t \langle w_i, \phi^n \rangle
$$

and

$$
(4.4) \quad \phi^{n+1} = \sum_{i=1}^{MM} \phi_i^{n+1} w_i + \sum_{MM+1}^{NS} \phi_{0i} w_i
$$

Exercise 1

Compute $b_{ij} = \int_{\Omega} w_i w_j \, dx$ and write a subroutine which computes $\phi^{n+1}$
by (4.3) (4.4), assuming $\phi^n$ known, by Choleski's factorization.

Then solve (4.1) (4.2) by using the subroutine N times. Check that if N is small the method diverges.

### Exercise 2

Solve (4.3), (4.4) and then (4.1) (4.2) by writing instead of

(4.3) (4.4)

$$ (4.5) \quad \phi_i^{n+1} = \phi_i^n - \Delta t (a(\phi^n, w_i) - b(w_i)) s(w_i) $$

$$ \text{where } s(w_i) = 3 \int_{\Omega} w_i \, dx $$

Formula (4.5) comes from a quadrature formula for $b_{ij}$. In some cases it is possible to show that this type of procedure leads to a proper scheme for (4.1) - (4.2). (see Raviart [1])

A fully *implicit* scheme would be (compare with (4.2))

$$ (4.6) \quad - \Delta \phi^{n+1} + \left(\frac{1}{\Delta t}\right) \phi^{n+1} = f^{n+1} + \frac{\phi^n}{\Delta t} $$

### Exercise 3

Solve (4.6) with FEMICG with $\Omega = [0, 1]^2, g = 0; \phi_o = 0$;

$$ f = 1 + \alpha\delta(t = \frac{1}{2}, x = (\frac{1}{2}, \frac{1}{2})) \quad \text{(Dirac function)} $$

$\alpha = 0 \text{ and } \alpha = 1$.

The *Crank-Nicholson* scheme would be

$$ - \frac{1}{2} \Delta \phi^{n+1} + \frac{1}{\Delta t} \phi^{n+1} = \frac{f^{n+1} + f^n}{2} + \frac{\phi^n}{\Delta t} + \frac{1}{2} \Delta \phi^n $$
### Exercise 4

Same as exercise 3 but with the Crank-Nicholson scheme, check that transitory states are bettwersimulated with the C-N scheme than with the implicit scheme. However the latter scheme is more stable.

For a good presentation of the time-derivative discretization problem see a Richmeyer & Morton [2]

## 2. The wave equation

$$ (4.7) \qquad \frac{\partial^2 \phi}{\partial t^2} - \Delta \phi = f \quad \text{in } \Omega \times ]0, T[, \quad T[ , \quad \phi|_{\Gamma} = 0 $$

$$ (4.8) \qquad \phi(0, x) = \frac{\partial \phi}{\partial t}(x, 0) = 0 $$

### Exercise 5

Use FEM1CG or FEM1CF to solve (4.7) (4.8) discretized by (4.18) chapter 4 of Lions [3] with $\lambda > 0$ (implicit)

Take $f = 1$, $\Omega = ]0, 1[^2$ and find the stability region by taking several values for w.

## 3. A non linear elliptic problem

The potential flow in a nozzle is governed by

$$ (4.9) \qquad \nabla \cdot \left[ \left( 1 - |\nabla \phi|^2 \right)^{2.5} \nabla \phi \right] = 0 \quad \text{in } \Omega \mid \Gamma_1 = \phi_0 \frac{\partial \phi}{\partial n} \mid \Gamma_2 = 0 $$
The speed of the gas is recovered
from $\phi$ by: $\mu = \nabla\phi$ in variational
form (4.9) is

<figure><img src="image_30.png" /><figcaption>Fig. 1</figcaption></figure>

$$
(4.10) \quad \int_{\Omega} (1 - |\nabla\phi|^2)^{2.5} \nabla\phi \nabla\psi \, dx = 0 \quad \forall \psi \in H_O^1(\Omega), \phi|_{\Gamma_1} = \phi_O
$$

<figure><img src="image_31.png" /></figure>

consider the problem

$$
(4.11) \quad \min_{\phi \in H_1^1(\Omega)} \{J(\phi) \mid \phi|_{\Gamma_1} = \phi_0\} \quad J(\Omega) = \int_{\Omega} (1 - |\nabla\phi|^2)^{3.5} dx
$$

Exercise 6

Show that J is strictly convex if $|\nabla\phi| < 1/4$

Show that the solution of (4.11) is a solution of (4.10).

(Show that $\lim_{\lambda \to 0} (J(\phi + \lambda\psi) - J(\phi)) = -7 \int_{\Omega} (1 - |\nabla\phi|^2)^{2.5} \nabla\phi \nabla\psi dx$ and change $\psi$ into $-\psi$). $\square$

<figure><img src="image_32.png" /></figure>

Therefore if $|\nabla\phi| < 1/4$ (subsonic flow) (4.10) has a unique
solution. One could use the SDM to solve (4.9). However it so
happens that the "quasi-linearisation" scheme below works very fast:

Compute $\phi^{n+1}$ from $\phi^n$ by solving

$$
(4.12) \quad \nabla \cdot \left[ (1 - |\nabla\phi^n|^2)^{2.5} \nabla\phi^{n+1} \right] = 0 \quad \phi^{n+1} \in H_O^1(\Omega), \phi^{n+1}|_{\Gamma_1} = \phi_O.
$$
### Exercise 7

Adapt FEM1CF or FEM1CG to solve (4.12) and solve (4.9) by doing 5 iterations of (4.12) with $\phi_0 = 0$ at the entrance of the nozzle and $\phi_0 = 0.8L$ at the exit. To discretize the nozzle use SQUARE and affine $x_2$-scaling. ($L$ is the length of the nozzle).

Check that if $\phi_0$ is increased so that $|\nabla\phi| > 1/4$, the method diverges.
For more details, see Periaux [4]. □

### 4. A linear elliptic constrained problem: Stokes flow

Stokes problem: $u(x) \in \mathbb{R}^n$, $p(x) \in \mathbb{R}$, $n = 2$ or $3$,

$$ (4.13) \qquad \Delta u = \nabla p + f, \nabla \cdot u = 0 \text{ in } \Omega, \quad u|_{\Gamma} = u_0 $$

has a variational formulation in $H^2(\Omega)$ (see Lions [3]). However a conformal discretization of $H^2(\Omega)$ by the FEM requires polynomial approximation of order 5 at least: this means the programming will be difficult. Instead we shall consider $\nabla \cdot u = 0$ as a constant. Then $p$ appears as the Lagrange multiplier associated with the constraint.
Therefore if we use the ALM we must solve

$$ (4.14) \qquad \min_{u \in U_h} \int_{\Omega} \left[ \frac{1}{2} |\nabla u^{n+1}|^2 - \rho^n \nabla \cdot u^{n+1} + f u^{n+1} \right] dx + \rho \int_{\Omega} |\nabla \cdot u^{n+1}|^2 dx $$

and take

$$ (4.15) \qquad p^{n+1} = p^n - \lambda \nabla \cdot u^{n+1} . $$
If the penalty method is used, one will solve

$$ (4.16) \quad \min_{u \in U_h} \int_{\Omega} [\frac{1}{2}|\nabla u|^2 + fu + \rho|\nabla \cdot u|^2]^2 dx $$

where $U_h = \{u = \sum_{i=1}^{MM} u_i w_i + \sum_{MM+1}^{NS} u_0 w_i\}$, if P1 elements are used.

### Exercise 8

Solve (4.16) with the conjugate gradient method and P1 elements; use some of the subroutines of Chapter 2. Take $\rho = 1000$ and $\Omega = [0, 1]^2$, $u_0 = 0$ except on the top of $\Omega$ where $u_0 = 1$ (driven cavity problem). For more details see Bercovier [5].

## 5. The incompressible stationary Navier Stokes equation

The problem is

$$ (4.17) \quad \nu \Delta u = \nabla \rho + u \nabla \cdot u = 0 \quad \text{in } \Omega \quad u|_{\Gamma} = u_0. $$

As is well-known, this equation has boundary layers when $u \to 0$. Any discretization scheme which has no nodes in the boundary layers will very likely not give a good result. Therefore to every discretization grid there will be a minimum allowable $\nu$. In fact, one can show that the real parameter is $R = \frac{u_0 L}{\nu}$ ($L$ = typical length scale of the flow). It turns out that it is very difficult to go beyond $R = 1000$.
For $R < 1000$ consider the algorithm: $\rho, \omega > 0$ large

$$ (4.18) \qquad \omega(u^{n+1} - u^n) + \nu\Delta\bar{u} = u^n\nabla u^n + \nabla p^n $$

$$ (4.19) \qquad \rho(p^{n+1} - p^n) = -\nabla \cdot u^{n+1} $$

with $\bar{u} = u^n$ (explicit) or $\bar{u} = u^{n+1}$ (implicit). This algorithm (uzawa's if $\bar{u} = u^n$) is a mixture of (4.14) - (4.15) and quasi-linearisation. For more details see Fortin [6] and Thomasset [7].

### Exercise 9

Solve (4.18) - (4.19) with $\bar{u} = u^n$, P1 finite elements for $u$ and $p$ (it is better if $u$ is P2, but it is more complicated) and quadrature approximations like (4.5). Take $\Omega = [0, 1]^2$ and $u_0 = 1$ on the top, zero elsewhere.

## 6. The Biharmonic problem

Consider the 4th order PDE

$$ (4.20) \qquad \Delta(\Delta\phi) = f \quad \text{in } \Omega $$

$$ (4.21) \qquad \phi|_{\Gamma} = \phi_0 \frac{\partial\phi}{\partial n}|_{\Gamma} = g. $$

This system admits a variational formulation in $H_0^2(\Omega)$ (see Lions [3]). However to discretize $H^2(\Omega)$ with conformal triangular elements, one needs polynomials of degree 5 at least. we present an alternative.

### 6.1 A splitting method

Note that if $\lambda = \Delta\phi|_{\Gamma}$ is known then (4.20) (4.21) split into two Dirichlet problems
$$ (4.22) \quad -\Delta\omega = f \quad \text{in } \Omega \qquad \omega|_{\Gamma} = \lambda $$

$$ (4.23) \quad -\Delta\phi = w \quad \text{in } \Omega \qquad \phi|_{\Gamma} = \phi_0 $$

Therefore let $\phi(\lambda)$ be the solution of (4.22) - (4.23) and let A be the mapping

$$ A : \lambda \to A\lambda = \frac{\partial\phi}{\partial n}(\lambda)|_{\Gamma} $$

Then the idea is to solve the linear system

$$ (4.24) \quad A\lambda = g. $$

One can show that A is affine symmetric coercive and that if (4.22) (4.23) are discretized by FEM with conformal element of degree k, the discrete analogue of A is positive definite. However (4.22) (4.23) (4.24) in discrete form is not equivalent to (4.20) (4.21) but to a mixed formulation of (4.20) - (4.21) and therefore k should be greater or equal to 2. For more details see Ciarlet Raviart [8] and Glowinski, Pironneau [9].
References for Chapter 4

[1] RAVIART, P. The use of numerical integration in FEM for solving Parabolic Equations in Topics on Numerical Analysis. J.H. Miller, ed. Academic Press 1972.

[2] RICHTMEYER, D. & MORTON, K. Difference methods for initial value problems. Interscience 1957.

[3] LIONS, J.L. These lecture notes.

[4] PERIAUX, J. 3-D analysis of compressible potential flows with the FEM. Int. J. for num. meth. in Eng. Vol. 9 1975.

[5] BERCOVIER, M. Thesis. University of Roven 1976.

[6] FORTIN, M. Thesis. University of Paris 6. 1972.

[7] THOMASSET. Application of FEM to the resolution of Navier Stokes equations. IRIA report No. 150, 1975.

[8] CIARLET, P.M., RAVIART, P.A. A mixed finite element method for the biharmonic equation. Proc. Symp. on Math. aspect of FEM. C de Boor ed. Academic Press N.Y. (1974).
# CHAPTER 5

## Methods for the Solution of some Variational Inequalities

In this section, we shall consider some examples of variational inequalities, some of which were studied by LIONS [1, Ch. 3]. (We shall restrict ourselves to the "symmetric" case).

### 5.1. The use of duality methods

Let $V = H_0^1(\Omega)$ where $\Omega$ is an open set of $\mathbb{R}^2$ and

$$J(v) = \frac{1}{2} a(v,v) - (f,v)$$

where

$$a(u,v) = \int_{\Omega} \nabla u \cdot \nabla v \, dx.$$

$$(f,v) = \int_{\Omega} f v \, dx \quad (f \in L^2(\Omega)).$$

We shall consider 2 variational inequalities for the 2 convex:

$$K_1 = \{ v \in H_0^1(\Omega) : v(x) \ge 0 \quad \text{a.e. } x \in \Omega \}$$

$$K_2 = \{ v \in H_0^1(\Omega) : |(\nabla v)(x)| < 1 \text{ are } x \in \Omega \}.$$

The corresponding variational inequalities are

$$(1) \quad a(u,v-u) - (f,v-u) \ge 0 \quad \forall v \in K_1$$

which constitutes the obstacle problem, and
$$ (2) \quad a(u, v-u) - (f, v-u) \ge 0 \quad \forall v \in K_2, $$

which constitutes the elasto-plastic torsion problem.

Note that only the convex K differs from one problem to the other. We shall see that the methods to be used are completely different.

We recall that (1), (2) are equivalent (since a is symmetric) to the following optimization problem:

$$ (3) \quad \begin{array}{ll} \text{Inf} & J(v) \\ v \in K & \end{array} $$

where K is equal to $K_1$ for the first problem, and to $K_2$ for the second one.

Of course, the first step for solving (3) consists of approximating $H_0^1(\Omega)$, which can conveniently be done by a finite element space, constituted of continuous piecewise affine functions over a triangulation $\mathcal{T}_h$ of $\Omega$ (see [1]). Let $V_h$ denote this finite element space, where $h$ denotes the diameter of the greatest triangle of $\mathcal{T}_h$.

We shall choose the approximate convex for both problems

$$ K_h^1 = K_1 \cap V_h $$

$$ K_h^2 = K_2 \cap V_h $$

The main difference for the solution of problems (1) and (2) is that it is easy to compute the projection over $K_h^1$, and not on $K_h^2$.
In fact, let $v_h \in V_h$, the projection of $v_h$ over $K_h^1$ is the function of $V_h$ which has the same value as $v_h$ at all the vertices where $v_h$ is positive and which takes the value zero at all the other vertices.

This means that we shall rather use the gradient algorithm with projection for solving the problem (1) (see §2.1) (or eventually the relaxation method, see LIONS [1, Ch. 3, §7.2]).

On the contrary, to solve the problem (2), we shall use either the penalty method (§2.2) or the duality method (§2.3, c)). As a matter of fact the experiments have shown the superiority of the duality method over the penalty method (see GLOWINSKI-LIONS-TREMOLIERES[2]) in the particular case of problem (2).

**Application of the gradient algorithm with projection to the solution of (1)**

First of all, we have to explain how to compute the gradient of $J$ considered as defined over $V_h$.

Introducing the stiffness matrix $A$ (see Ch. 3), which is a $N \times N$ matrix, where $N$ is the dimension of $V_h$ (i.e., the number of interval vertices of the triangulation $T_h$), we know that

$$J(v_h) = \frac{1}{2} V_h^T A V_h - V_h^T b$$

where $V_h$ is the column vector of the values of the function $v_h$ at the internal vertices of $T_h$, and $b$ the right hand side vector defined in Ch.3, formula (6).
We see that if we choose as scalar product of $V_h$ the Euchidian scalar product of $R^N$, then the gradient of $J$
will be $J'(v_h) = AV - b$
(see exercise 2, Ch. 1) (Note that this define an element of $V_h$).

The algorithm will then be the following:

. $v^0 \in \mathbb{R}^N \text{ and } \rho > 0 \text{ being given;}$

(4)

. $v^n \in \mathbb{R}^N \text{ being given by induction;}$

. compute $G^n = AV^n - b$

. set $v_i^{n+1} = \max(0, v_i^n - \rho G_i^n)$ for $i = 1, \dots, N$. □

**Remark:** Note that if $K$ is equal to the whole space $H_0^1(\Omega)$, the algorithm (4) is a simple gradient algorithm which is known to be slow if the condition number of the matrix $A$ was big (see Ch. 1, §3), which is the case when $h$ is very small (N very big).

In such a case it is convenient, either to adapt the conjugate gradient algorithm to this case with constraints, (cf. TREMOLIERES [3]), or to use Uzawa's algorithm (see Ch. 1, §2.3). For the last one we have to solve a linear system of equations at each step, the matrix of which is $A$. This can be conveniently done by Gaussian elimination. □

**Application of the duality method to problem (2)**

As $|\text{grad } v_h|$ is piece wise constant, we have as many constraints
as triangles in $\bar{T}_h$. In a very natural way, we first define for
$\mu \in L^\infty(\Omega)$ the Lagrangian

$$L(v,u) = J(v) + \int_{\Omega} \mu(|\nabla v|^2 - 1) \, dx.$$

We check that if $\mu_h$ is a piecewise constant function and $v_h \in V_h$, then $L(v_h, \mu_h)$ is the equivalent of the Lagrangian introduced for inequality constraints in Ch. 1, §2.3.

We can then use the Uzawa algorithm.

Note that to minimize $L(v,\mu)$ in $v$, $\mu$ being fixed, we just have to minimize the quadratic functional

$$\frac{1}{2} \int_{\Omega} (1 + \mu(x)) |\nabla v|^2 dx - \int_{\Omega} f v dx.$$

Then, $\mu_h$ being fixed (in the Uzawa algorithm) we have to solve a linear system of equations

$$A(\mu_h) V = b$$

where $A(\mu_h)$ is the stiffness matrix associated to the bilinear form

$$\int_{\Omega} (1 + \mu_h) \nabla u_h \cdot \nabla v_h \, dx.$$

(Note that the stiffness matrix $A(\mu_h)$ is the result of an assembling of all the standard stiffness submatrices $A^K$ multiplied by a scale factor $(1 + \mu_h^K)$ if $\mu_h^K$ is the (constant) value of $\mu_h$ on the element K).
Uzawa's algorithm can be written in this special case as
follows:

$p > 0$ and $\mu^0$ being given (piecewise constant);

• $\mu^n$ being given (piecewise constant) by induction.

• compute the stiffness matrix $A(\mu^n)$ associated to $\mu^n$;

• solve the linear system

$$A(\mu^n) v^n = b.$$

• compute $\mu^{n+1}$: on each element K, the restriction of $\mu^{n+1}$ is given by

$$\mu|_{K}^{n+1} = \max (0, \mu|_{K}^{n} + \phi (|\nabla v|^{n} - 1))$$

where $v^n$ is the function of $v_h$ associated to the column vector $v^n$. □

## 5.2. Solution of some problems related to plasticity.

### a) Linear elasticity

We shall first consider the case of linear elasticity.

Consider a continuous medium $\Omega \subset \mathbb{R}^n$ fixed over a part $\Gamma_u$ of its boundary, and submitted to both pressure actions $F_i$ at the other part $\Gamma_F$ of the boundary, and to body forces $f_i$ inside $\Omega$. The stress tensor $\sigma$ (which is symmetric) must satisfy the equilibrium equation
$$
(5) \quad \sum_{j=1}^{n} \frac{\partial \sigma_{ij}}{\partial x_j} + f_i = 0 \quad \text{in } \Omega, \quad 1 \le i \le n
$$

$$
(6) \quad \sum_{j=1}^{n} \sigma_{ij} n_j = F_i \quad \text{in } \Gamma_F \qquad 1 \le i \le n
$$

where $n_j$ denote the components of the unit normal vector to $\Gamma_j$ outwards directed.

Moreover, the displacement field must satisfy the boundary
condition.

$$
(7) \quad u_i = 0 \quad \text{on } \Gamma_u, \quad 1 \le i \le n.
$$

We shall look for the displacements in the Sobolev Space

$$
V = \{ v \in (H^1(\Omega))^n \text{ such that } v = 0 \text{ on } \Gamma_u \}
$$

and the stress field in

$$
H = \{ \tau : \tau_{ij} \in L^2(\Omega), \tau_{ij} = \tau_{ji} \text{ for } i, j = 1, \dots, n \}
$$

which we provide with the scalar product

$$
(\sigma, \tau) = \int_{\Omega} \sum_{i,j=1}^{n} \sigma_{ij} \tau_{ij} \, dx
$$

associated to the norm $|| \cdot ||$.

<footnote>(*) n = 2 or 3 in the applications.</footnote>
Let us introduce the partial derivative operator

$$ \varepsilon : V \to H $$

defined by

$$ \varepsilon_{ij}(v) = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \quad 1 \le i,j \le n. $$

The tensor $\varepsilon(u)$ associated to the displacement field $u$ is called strain tensor.

In the case of linear elasticity, there is a linear relation between $\sigma$ and $\varepsilon(u)$, of the following type.

$$ (8) \quad \sigma_{ij} = \lambda \sum_{k=1}^{n} \varepsilon_{kk}(u) \delta_{ij} + 2\mu \varepsilon_{ij}(u) $$

where $\delta_{ij}$ denotes the Kronecker tensor (*) and $\lambda$ and $\mu$ are the Lamé coefficients.

For the sake of simplicity, we shall assume that $\lambda = 0$ and $\mu = \frac{1}{2}$ (this does not change the fundamental nature of the problem), so that

$$ (9) \quad \sigma = \varepsilon(u) $$

In order to find a variational formulation of the problem (5) (6)(7)(9), we prove

<footnote>(*) $\delta_{ij} = 1$ if $i = j$ and 0 otherwise</footnote>
**Lemma 1 (Virtual work theorem):** Let L be the linear form over V:

$$L(v) = \sum_{i=1}^{n} \left( \int_{\Omega} f_i v_i dx + \int_{\Gamma_F} F_i v_i d\gamma \right)$$

then if $\sigma$ satisfies (5)(6), then one has

$$(10) \quad (\sigma, \varepsilon(v)) = L(v) \quad \forall v \in V.$$

**Proof:** We use Green's formula (and the symmetry of $\sigma$)

$$\int_{\Omega} \sigma_{ij} \varepsilon_{ij}(v) \, dx = \int_{\Gamma_F} \sigma_{ij} n_j v_i \, d\sigma - \int_{\Omega} v_i \frac{\partial \sigma_{ij}}{\partial x_j} \, dx. \quad \square$$

The linear elasticity problem is then reduced to the variational problem (replace (9) in (10)):

Find $u \in V$ such that

$$(11) \quad (\varepsilon(u), \varepsilon(v)) = L(v) \quad \forall v \in V.$$

(To obtain the more general problem, one would have to replace (8) in (10)).

Note that to prove the existence of a solution to (11), we have to prove first that the bilinear form $(\varepsilon(u), \varepsilon(v))$ is elliptic. This is a non trivial proof which is performed in DUVAUT-LIONS [4]

b) Solution of a linear elastic problem

Assume that the problem is 2-dimensional (n = 2; this is the case either for a plane stress or a plane strain problem).
We still introduce a triangulation of $\Omega$, and define $V_h$ as the subspace of $V$ constituted with vector valued functions each component of which is continuous and piecewise affine.

The dimension of the space $V_h$ is here 2 times the number of vertices which are not on $\Gamma_u$.

The stiffness submatrix of an element will be now 6 x 6 instead of 3 x 3, since there are 2 components per vertex.

Apart from this, the computations are handled in exactly the same manner as in chapter 3 for a Dirichlet problem.

### c) Elasto-plasticity

Let us consider first the case of the Hencky law:

(see DUVAUT-LIONS [4]) (which is to replace (8) or (9)).

$$ (12) \quad (\sigma, \tau - \sigma) - (\varepsilon(u), \tau - \sigma) \ge 0 \qquad \forall \tau \in K $$

where $K$ is the plasticity convex, namely

$$ K = \{ \tau \in H : F(\tau) \le 0 \quad \text{a.e. on } \Omega \}. $$

and $F$ the plasticity criterion.

Let us introduce the set of statically admissible stress fields (i.e. satisfying (5) and (6)):

$$ E = \{ \tau \in H : (\tau, \varepsilon(v)) = L(v), \quad \forall v \in V \}. $$
We see that for $\tau \in K \cap E$, (12) reduces to

$$ (\sigma, \tau - \sigma) \ge 0, $$

in such a way that $\sigma$ is the solution (if any) of the
optimization problem

$$ (13) \quad \inf_{\tau \in K \cap E} \frac{1}{2} \|\tau\|^2. $$

If $K \cap E$ is not empty, (i.e., if the external loads $f$ and $F$ are not too big), then there exists a unique solution $\sigma$ to problem (13).

However it is more difficult to define a displacement field
associated to $\sigma$ (see MERCIER [5])

The approximation of problem (13) is done in the following
way. Let $V_h$ be defined as for the elastic problem (cf
section b)), and $H_h$ be the subspace of $H$ of piecewise
constant tensor fields, we choose

$$ E_h = \{\tau \in H_h : (\tau, \varepsilon(v)) = L(v), \quad \forall v \in V_h\} $$

and approximate (13) by

$$ (14) \quad \inf_{\tau \in K \cap E_h} \frac{1}{2} \|\tau\|^2 $$

Let $\sigma_h$ denote the solution of (14), we have the convergence

result: $\sigma_h \to \sigma$ as $h \to 0$.

(see MERCIER [5])
To solve problem (14), we notice that this is very similar to problem (26) chapter 1.

We then use the gradient method with auxiliary operator to solve the dual problem, which is in this case equivalent to find $u_h$ solution of

$$ (P_K(\varepsilon(u_h)), \varepsilon(v)) = L(v) \quad \forall v \in V_h, $$

where $P_K$ denotes the projection operator onto $K$.

The application of the gradient method with auxiliary operator to this problem gives the algorithm (see [5] for details):

1. $u^0 \in V_h$ and $\rho > 0$ being given;

2. $u^n \in V_h$ being given, by induction;

3. compute $u^{n+1}$ solution of the (linear elastic) problem

$$ (\varepsilon(u^{n+1}), \varepsilon(v)) = (\varepsilon(u^n), \varepsilon(v)) - \rho[(P_K(\varepsilon(u^n), \varepsilon(v)) - L(v))] $$

$$ \forall v \in V_h. \quad \square $$

This algorithm is shown to converge for $0 < \rho < 2$ by using Theorem 2, Chapter 1.

### d) Quasi-static elasto-plastic problem

To reproduce unloading phenomena, it makes more sense to use the Prandtl-Reuss law.

This law concerns the case of an external load slowly varying
with time. (The linear form L introduced at Lemma 1 is then
depending on time).

We denote by $\dot{u}$, and $\dot{\sigma}$ the partial derivatives with
respect to $t$ of $u(t,x)$ and $\sigma(t,x)$.

The Prandtl-Reuss law can be written as follows (see [4])

$$
(15) \quad (\dot{\sigma}, \tau - \sigma) - (\varepsilon(\dot{u}), \tau - \sigma) \ge 0 \quad \forall \tau \in K, t \ge 0;
$$

this relation is completed by the one expressing
that $\sigma(t)$ should be in equilibrium for each $t \ge 0$:

$$
(\sigma(t), \varepsilon(v)) = \langle L(t), v \rangle \quad \forall v \in V
$$

(We denote by $\langle \cdot, \cdot \rangle$ the duality pairing between $V$ and $V'$).

We complete with the initial conditions:

$$
\sigma(0) = \sigma^0; \quad u(0) = u^0.
$$

Let

$$
E(t) = \{ \tau \in H : (\tau, \varepsilon(v)) = \langle L(t), v \rangle \quad \forall v \in V \}
$$

in a similar way to the stationary case, we can eliminate
$\dot{u}$ in (15) by choosing $\tau \in E(t)$. We get that $\sigma$ is solution
of the following problem

$$
(18) \quad (\dot{\sigma}, \tau - \sigma) \ge 0 \quad \forall \tau \in K \cap E(t)
$$
The existence of a solution to problem (18) has been proved by JOHNSON [6], when $K \cap E(t)$ is non empty, and that $L(t)$ is not too much varying.

To solve the problem (18), we use an implicit scheme:

$$ (19) \quad (\sigma^{n+1} - \sigma^n, \tau - \sigma^{n+1}) > 0 \quad \forall \tau \in K \cap E(t^{n+1}) $$

which is in fact a problem very close to the stationary problem (13), and that can be solved by the same methods.

Note that the time step $\Delta t$ has disappeared in (19), and this is due to the fact that a quasi-static evolution is supposed to be infinitely slow.
References for Chapter 5

[1] LIONS, J.L.: These lecture notes.

[2] GLOWINSKI, R., LIONS, J.L., TREMOLIERES, R.

[3] TREMOLIERES, R.: Inequations variationnelles: existence, approximation, resolution. These, Paris, 1972.

[4] DUVAUT, G., LIONS, J.L., Lesinequations en niecauique et en physique, Dunod, Paris, 1972.

[5] MERCIER, B.: Sur la Theorie et l'analyse numerique de problinues de plasticite, These, Paris, 1977.

[6] JOHNSON, C.: "Existence theorems for plasticity problems", J. Maths pures et Appliquees 55 (1976) 431-444.
## Chapter 6

### Optimal Control and Optimum Design

#### O. Preliminaries

We begin with a brief introduction to numerical optimization in Hilbert spaces.

Let $H$ be a Hilbert space with $\langle \cdot, \cdot \rangle$ its scalar product and $|| \cdot ||$ its norm.

#### Example

$$ H = L^2(\Omega) \quad \langle \phi, \psi \rangle = \int_{\Omega} \phi(x) \psi(x) dx $$

$$ H = H^1(\Omega) \quad \langle \phi, \psi \rangle = \int_{\Omega} [\nabla\phi(x)\nabla\psi(x) + \phi(x)\psi(x)] dx $$

Given $J: H \to R$, we shall consider the problem

$$ \min J(u) \\ u \in H $$

#### Definition

1. $J$ is F-differentiable at $u$ if there exists $J: H \to H'$ (dual of $H$) such that

$$ \forall h \in H \quad \lim_{\lambda \to 0} \frac{1}{\lambda} (J(u + \lambda h) - J(u)) = (J'(u), h) $$

where $(\cdot, \cdot)$ is the product of duality between $H$ and $H'$.
2. If $J$ is F-differentiable, the gradient of $J: \mathrm{grad} \, J: H \to H$ is defined by

$$ \langle \mathrm{grad} \, J(u), h \rangle = (J'(u), h) $$

i.e. $\mathrm{grad} \, J(u)$ is the unique element of $H$ such that

$$ J(u + \lambda h) - J(u) = \lambda \langle \mathrm{grad} \, J(u), h \rangle + O(\lambda) $$

$$ \text{with } \lim_{\lambda \to 0} \frac{1}{\lambda} O(\lambda) = 0. $$

Now we have a formula identical to (3) Ch. 1, therefore we can write the method of steepest descent or any other method of Chapter 1 to solve our problem. For example:

**Algorithm:** SDM

Step 0: Choose $u^0 = 0$, set $i = 0$;

Step 1: Set $d^i = - \mathrm{grad} \, J(u^i)$;

Step 2: Compute $\lambda^i \in R$ solution of
$$ \min J(u^i + \lambda^i d^i); $$
$\lambda$

Step 3: Set $u^{i+1} = u^i + \lambda^i d^i$ and go back to 1. □

# 1. Resolution of an Optimal Control Problem

To illustrate the technique we begin with a very simple example.

## 1.1. Example

What is the electrostatic charge distribution in $\Omega$ which will give a desired potential distribution $\phi_d$ in $\Omega$; i.e. find $f$ such that
the solution $\phi$ of

$$ (6.1) \quad \Delta\phi = f \text{ in } \Omega \quad \phi|_{\Gamma} = 0 $$

is equal to $\phi_d$ in $\Omega$. One way of solving

the problem is to find $\tilde{f}$ solution of the problem

$$ (6.2) \quad \begin{array}{ll} \min & J(f) \\ \text{subject to} & f \in L^2(\Omega) \end{array} $$

where

$$ (6.3) \quad J(f) = \int_{\Omega} |\nabla(\phi(f) - \phi_d)|^2 dx $$

with $\phi(f)$ solution of (6.1). Problem (6.2) is an optimal control problem for a distributed parameter system.

## 1.2. Solution of (6.2) by SDM

First we must find grad $J(f)$. By definition, since the scalar product of $L^2(\Omega)$ is

$$ (6.4) \quad \langle f, g \rangle_{L^2(\Omega)} = \int_{\Omega} f(x)g(x)dx. $$

$J$ is the function of $L^2(\Omega)$ such that

$$ (6.5) \quad J(f + \delta f) - J(f) = \int_{\Omega} \operatorname{grad} J(f)(x) \cdot f(x) dx + \int_{\Omega} \delta f (\| \delta f \|) \wedge \delta f \in L^2(\Omega) $$

where $o(\cdot)$ is such that

$$ (6.6) \quad \lim_{\lambda \to 0} o(\lambda)/\lambda = 0 $$
Now from (6.3)

$$ (6.7) \quad J(f + \delta f) - J(f) = \int_{\Omega} [ |\nabla(\phi(f + \delta f) - \phi_d)|^2 - |\nabla(\phi(f) - \phi_d)|^2 ] dx . $$

Let

$$ (6.8) \quad \delta\phi(x) = \phi(f + \delta f)(x) - \phi(f)(x). $$

Then

$$ (6.9) \quad J(f + \delta f) - J(f) = 2 \int_{\Omega} \nabla(\phi(f) - \phi_d) \cdot \nabla \delta\phi dx + \int_{\Omega} |\nabla \delta\phi|^2 dx. $$

It is shown in Lions [2] that $\int_{\Omega} |\nabla \delta\phi|^2 dx = o(\| \delta f \|)$. By definition

(i.e. (6.1)) $\delta\phi$ is the solution of

$$ (6.10) \quad \Delta\delta\phi = \delta f \text{ in } \Omega \quad \delta\phi|_{\Gamma} = 0. $$

In variational form (6.10) is

$$ (6.11) \quad \int_{\Omega} \nabla \delta\phi \nabla v dx = \int_{\Omega} \delta f v dx \quad \forall v \in H_0^1(\Omega). $$

Let $v = \phi(f) - \phi_d$ in (6.11), then

$$ (6.12) \quad \int_{\Omega} \nabla \delta\phi \nabla (\phi(f) - \phi_d) dx = \int_{\Omega} \delta f (\phi(f) - \phi_d) dx $$

Therefore (6.9) becomes
$$ (6.13) \quad J(f + \delta f) - J(f) = 2 \int_{\Omega} \delta f (\phi(f) - \phi_d) dx + o(||\delta f||) $$

and by identification with (6.5) we find that

$$ (6.14) \quad \mathrm{grad} \, J(f)(x) = 2(\phi(f)(x) - \phi_d(x)). $$

### Remark

In general it is not possible to proceed from (6.12) to (6.13). For example if J was

$$ (6.15) \quad J(f) = \int_{\Omega} |\phi(f) - \phi_d|^2 dx . $$

Then (6.9) would have been

$$ (6.16) \quad J(f + \delta f) - J(f) = 2 \int_{\Omega} (\phi(f) - \phi_d) \delta\phi dx + o(||\delta f||) . $$

One must introduce an adjoint state $\psi(x) \in H_0^1(\Omega)$ solution of

$$ (6.17) \quad \int_{\Omega} \nabla\psi\nabla v \, dx = 2 \int_{\Omega} (\phi(f) - \phi_d) v \, dx \quad \forall v \in H_0^1(\Omega) $$

$$ \text{i.e.} \quad -\Delta\psi = 2(\phi - \phi_d), \quad \psi|_{\Gamma} = 0. $$

Then by letting $v = \delta\phi$ in (6.17) and $v = \psi$ in (6.11) one finds that

$$ (6.18) \quad J(f + \delta f) - J(f) = \int_{\Omega} \psi(x) \delta f(x) dx \quad \forall \delta f \in L^2(\Omega) $$

therefore $\mathrm{grad} \, J(f)(x) = \psi(x)$
We can write the SDM algorithm for the solution of (0.2).

### Algorithm 1

0. choose $f^0 \in L^2(\Omega)$ set $i = 0$

1. compute $\phi^i$ solution of $-\Delta\phi^i = f^i$, $\phi^i|_\Gamma = 0$

2. set $f^i(\lambda) = f^i - 2\lambda(\phi^i - \phi_d)$, an compute $\lambda^i$ solution of
$$ \min_{\lambda} \int_{\Omega} |\nabla|\phi^i(f^i(\lambda)) - \phi_d|^2 dx $$

3. set $f^{i+1} = f^i(\lambda^i)$, $i = \ell+1$ and go back to 1.

### 4. Implementation of Algorithm 1

One possibility would be to discretize by FEM the equations of algorithm 1. In all likelihood this would yield on very poor algorithm because it would diverge when the discretization error on $\phi^i$ is of the order of $\phi^i - \phi_d$. It is much better to re-derive everything for the discrete optimal control problem. Therefore let us solve

$$ (0.20) \quad J_h(f_h) = \int_{\Omega} |\nabla(\phi_h(f_h) - \phi_{dh})|^2 dx $$

and $\phi_h(L_h)$ is the solution of

$$ (0.21) \quad \int_{\Omega} \nabla\phi_h \cdot \nabla v_h = \int_{\Omega} f_h v_h dx \quad \forall v_h \in V_{oh} $$

where, for example, $\tilde{L}_h = \{L_h | f_h|_T = \text{etc.} \quad \forall T \in T\}$
and $V_{oh} = \{v_h \in C^0(\Omega) \mid v_h|_T \text{ affine } \forall T \in T, v_h|_T = 0\}$.

If $\phi_{dh} \in V_{oh}$ then the computation of $\text{Grad} J_h(f_h)$ is done exactly as in the continuous case. We have

$$ (6.22) \quad J_h(f_h + \delta f_h) - J_h(f_h) = 2 \int_{\Omega} \delta f_h (\phi_h(f_h) - \phi_{dh}) dx + o(\| \delta f_h \|). $$

Therefore on each triangle T

$$ \delta f_h|_T = 2 \int_T [\phi_h(f_h) - \phi_{dh}]dx/mes T \triangleq \overline{\phi_h - \phi_{dh}}. $$

by definition of the symbol $\triangleq$.

### Algorithm 2

0. choose $f_h^0 \in L_h$ (for example $f_h^0 = 0$) for all $i = 0$

1. compute $\phi_h^i$ by (6.21)

2. compute (by the DPF for example, although in our case $\lambda^i$ can be computed explicitly) an approximation $\overline{\lambda^i}$ of the solution of

$$ (6.23) \quad \min_{\lambda} \int_{\Omega} |\nabla[\phi_h(f_h^i - 2\lambda(\overline{\phi_h^i - \phi_{dh}}))(x) - \phi_{dh}(x)]|^2 dx $$

3. Det $f_h^{i+1} = f_h^i - 2\lambda^i (\overline{\phi_h^i - \phi_{dh}})$, $i = i + 1$ and go back to 1.

### Remark

Note that the *clearest step is step 2* because each trial $\lambda$ in DPF cost the computation of $\phi_h(\lambda)$, i.e. the resolution of a PDF.

### Programming of Algorithm 2

We use the subroutine FEMICG for step 1 and DPF for step 2.
RELEASE 2.0

MAIN

DATE = 77067

15/26/37

C PROGRAMMING INITIALIZATION

DIMENSION Q1(2,3)
COMMON/A/UD(25),U1(25),U(25),F(32),F1(32),H(25),Q(2,25),G(1),DJ(25)
1)
COMMON/B/NS,NT,ME(3,32),NG(25)
CALL SQUARE(5.5,NS,Q,NT,ME,NG,MUD,NF,2)
R00=0.
R0=.001
DO 4 I=1,NS
UD(I)=0.
4 U(I)=1.

C
C STEP 0

C
DO 1 J=1,NT
1 F(J)=1.
I=0

C
C STEP 1

3 CONTINUE
CALL FEM1CG( 1.0..F,G,1.0,Q,ME,NG,NS,NT,U,50,.1E-9,1,H,DJ)
PRINT 101,(F(L),L=1,32)
PRINT 101,(U(L),L=1,25)
101 FORMAT(5(3X,E9.2))

C
C STEP 2

C
C COMPUTES J(F) BY 6.3
CALL CPS(R00,C)
CALL DPF(R0,C)

C
C STEP 3

C
C COMPUTES ALSO GRAD2=LIMIT(J(F(R0))-J(F(0)))/R0
GRAD2=0.
DO 2 J=1,NT
UUD=0.
DO 6 K=1,3
MEK=ME(K,J)
DO 7 L=1,2
7 Q1(L,K)=Q(L,MEK)
6 UUD=(U(MEK)-UD(MEK))/3. +UUD
AIR=.5*(Q1(1,1)-Q1(1,3))*(Q1(2,2)-Q1(2,3))-.5*(Q1(1,2)-Q1(1,3))*(Q1
11(2,1)-Q1(2,3))
IF(AIR>.0.)AIR=-AIR
GRAD2=GRAD2+4.*UUD*UUD*AIR
2 F(J)=F(J)-2.*RO*UUD
I=I+1

C
C STOP RULE

C
PRINT 100.C,GRAD2,RO
100 FORMAT(1X,'CRIT= ',E13.6,'GRAD2=' ,E9.2,' R00PT= ',E9.2)
IF(I,LE,10,AND,C,GT,.0001)GO TO 3
STOP
END

SUBROUTINE CPS(R0,C)

C COMPUTES SUM ON ΩΜΕGA OF GRAD(U(F-2RO(U-UD)BAR))-UD) **2
COMMON/A/UD(25),U1(25),U(25),F(32),F1(32),H(25),Q(2,25),G(1),DJ(25
1)
COMMON/B/NS,NT,ME(3,32),NG(25)
DO 5 J=1,NS
5 U1(J)=U(J)
DO 2 J=1,NT
UUD=0.
DO 4 K=1,3
MEK=ME(K,J)
4 UUD=UUD -2.*RO*(U(MEK)-UD(MEK))/3.
2 F1(J)=F(J)+UUD
CALL FEM1CG( 1.0..F1,G,1.0,Q,ME,NG,NS,NT,U1,20,.1E-7,1,H,DJ)
DO 3 I=1,NS
3 U1(I)=U1(I)-UD(I)
CALL AHH(Q,ME,NS,NT,NG, 1.0..U1,C)
RETURN
END
**Exercise 1:**

Choose a $\phi_d$ and run this program by completing it with

- the Dimension statement

- the Construction of a triangulation : use CERCLE

- the Printing statements for the results. □

**Exercise 2:**

Same as above but now $J_h$ is given by (6.15). □

**Exercise 3 (1h)**

Same as above but now $L_h$ is replaced by $V_h$

$$V_h = \{v_h \in C^0(\Omega) \mid v_h|_T \text{ affine } \forall T \in T\}. \quad \square$$

**Exercise 4 (10h)**

Same as above but (6.21) is integrated with $Q_1$ elements. □

## 2. Some Examples of Applications of Optimal Control

### 2.1. Control of the temperature in an oven

What is the temperature distribution $\tau_o(t)$ ($t$ is the time) of $\Gamma$ that will bring a piece of metal $\Omega$, $\partial\Omega = \Gamma$, to its melting temperature $\tau_m$ at time $T$? One may solve
$$
(6.25) \quad \min_{\tau_o \in Q} \left\{ \int_0^T \left( \int_\Omega |\tau - \tau_m|^2 dx dt \right) \left| \frac{\partial \tau}{\partial t} - \Delta \tau = 0 \text{ in } \Omega, \tau(t)|_\Gamma = \tau_o(t) \right. \right\} \\
\qquad \left. \tau(0, x) = 0 \quad x \in \Omega \right\}
$$

Exercise 5

Solve (6.25) by taking $Q = \{\tau_o\}$ piecewise constant in $t$ on $N$
intervals of $[0, T]$ of size $\Delta t = T(N+1)$; piecewise linear and continuous
on $\Gamma$. Discretize $\frac{\partial \tau}{\partial t} - \Delta \tau$ by the FEM with $P1$ elements and by an
implicit scheme in time. Derive the optimally conditions (see Lions [1])
and use the SDM. (See also YVON et al [2]). $\square$

2.2. Integration of non-linear PDF

The small disturbance equation for potential compressible
flow (see Landau-Lifshitz [3]) is

$$
\Delta\phi = \frac{\partial}{\partial x} \left[ \left( \frac{\partial\phi}{\partial x} \right)^2 \right] + \frac{\partial^2\phi}{\partial x^2} \text{ on } \Omega \quad \left. \frac{\partial\phi}{\partial n} \right|_{\Gamma_1} = g \quad \phi|_{\Gamma_2} = \phi_o \quad \Gamma_1 \cup \Gamma_2 = \Gamma
$$

This equation is elliptic if $\frac{\partial \phi}{\partial x} < 0$ and hyperbolic if $\frac{\partial \phi}{\partial x} > 0$. The
quasilinearisation method fails in this case, however one can solve
the optimal control problem

$$
(6.26) \quad \min_{u \in U} \left\{ \int_{\Omega} \left( u - \frac{\partial \phi}{\partial x} \right)^2 dx \left| \begin{array}{l} \Delta \phi = \frac{\partial}{\partial x} (u^2 + u) \frac{\partial \phi}{\partial n} \Big|_{\Gamma_1} = g \quad \phi|_{\Gamma_2} = \phi_o \end{array} \right. \right\}
$$
### Exercise 6

Solve (6.26) with $\Omega = \frac{1}{2}$ nozzle, $\phi$ integrated by the FEM with P1 elements (the right member should be integrated by part) and $U = L_h$. Use the SDM but compute explicitly the optimum step size $\lambda$. (See Glowinski et al [3]). □

### Exercise 7 (30h)

Same as above but with the Navier Stokes equation where the equation $\nabla \cdot u = 0$ is penalized. □

### 3. Optimum Design

It is beyond the scope of this lecture note to solve completely an optimum design problem. However, for research student who would like to try we shall give some clues:

Consider the problem

$$ (6.27) \quad \min_{\Omega \in \Theta} \left\{ \int_E \left| \phi - \phi_d \right|^2 dx \left| \begin{array}{l} -\Delta \phi + \phi = f \quad \text{in } \Omega \\ \partial_n \phi = 0 \end{array} \right. \right\} $$

and assume for simplicity that

$$ U = \{ \Omega \text{ open } \bar{E} \subset \Omega \} $$

$$ f = 0 \text{ in } \Omega - E. $$

The optimality conditions can be found by taking $\Omega'$ close to $\Omega$ with $\Omega' = \{x + \alpha(x)n(x) | x \in \Omega, n(x) = \text{normal of } \partial\Omega \text{ at } x\}, \alpha \text{ small.}$

One finds that if
$$
(6.28) \quad J(\Omega) = \int_E |\phi - \phi_d|^2 dx
$$

$$
(6.29) \quad J(\Omega') - J(\Omega) = - \int_{\partial\Omega} \nabla\phi \cdot \nabla \text{pad}\Gamma + o(|\alpha|)
$$

where p is the solution of

$$
(6.30) \quad \int_{\Omega} (\nabla p \cdot \nabla v + pv) \, dx = 2 \int_{E} (\phi - \phi_d) v \, dx \quad \forall v \in H^1(\Omega)
$$

$$
\begin{flalign*}
& \text{(this is because, where } \Omega' > \Omega: && \int_{\Omega} (\nabla \delta\phi \cdot \nabla v + \delta\phi v) \, dx + \\
& && \int_{\Omega' - \Omega} (\nabla \phi' \cdot \nabla v + \phi' v) \, dx = \int_{\Omega' - \Omega} f v \, dx \quad \forall v \in H'(\Omega') &&
\end{flalign*}
$$

To obtain the discrete analogue of (6.29), one should
proceed as follows:

(1) show that if $w^i$ the affine function on T
such that $w^i(x^i) = 1$, $w^i(x^k) = 0$, $w^i(x^j) = 0$
and if $x^k$ moves by $\delta x^k$ then $w^i$ becomes
$w^i + \delta w^i$ with

$$
(6.31) \quad \delta \mathbf{w}^i(x) = - \mathbf{w}^k(x) \xrightarrow{\nabla \mathbf{w}^i} \delta \mathbf{x}^k \quad \forall x \in T
$$

(use (2.26) and (2.27)).

(2) Then is $\phi(\Omega) = \sum \phi^i w^i$ and $\phi(\Omega') = \sum \phi'^i w'^i$ by definition, if $\phi$ is integrated with P1 elements

<figure><img src="image_33.png" /><figcaption>Fig. 1</figcaption></figure>
(4) By introducing the adjoint state p as in (6.30), show that

$$
(6.34) \quad J(\Omega') - J(\Omega) = \sum_k \int_{\Omega} (\nabla\phi \cdot \alpha_k \nabla p \cdot \nabla w^k + \nabla\phi \cdot \nabla w^k \nabla p \cdot \alpha_k)dx \\
\qquad - \sum_k \sum_{j \in J_k} \sum_{i \in I_{jk}} \frac{1}{2} \int_{L_i} (\nabla\phi \nabla p) (\alpha_k^T \vec{\hbar}) d\Gamma .
$$

The last integral is the discrete analogue of (6.29) and the first
two come from the discretization, they tend to zero when the
discretization parameter h tends to zero because $\frac{\partial\phi}{\partial n}$ and $\frac{\partial p}{\partial n}$ are zero
on Γ. Show by Green's formula that the last integral is equal to
$$ - \int_{\Omega} \nabla\phi \cdot \nabla p \, \nabla w_k \cdot \alpha_k dx . $$

Now the interested reader can program the SDM. Note that
great care should be taken so that the triangulations move in a
continuous fashion. The best way to insure it is to assign each node
to stay on prescribed lines, see Fig. 2. For more details, see
Marrocco-Pironneau [4]. When Γ moves in θ, one way is to keep the
nodes on the vertical lines and to move the interior nodes of a
quantity proportional to the displacement of the corresponding node
of Γ on the vertical line.

<figure><img src="image_34.png" /><figcaption>Fig. 2</figcaption></figure>
References for Chapter 6

[1] LIONS, J.L.: Optimal control of distributed parameter systems.
Springer 1969.

[2] YVON, J.P.: These d'etat universite de Paris 6, 1972.

[3] GLOWINSKI, R. et PIRONNEAU, O.: Calcul d'econlement transsenique
en elements Liris par le controle optimal. Proc. IRIA conf.,
Dec. 1975.

[4] MARROCCO, A. et PIRONNEAU, O.: To appear.