Implementación del Método de los Elementos Finitos

Por fijar ideas y para facilitar las explicaciones, vamos a describir el proceso de resolución de un problema por el método de los elementos finitos (MEF) sobre un problema concreto.

$\Omega\subset\mathbb{R}^2$ el abierto acotado representado en la figura.

Consideramos el siguiente problema:

\begin{matrix} (P) & {\begin{cases} - Δ u + u = f & en Ω, \\ u = u_{0} & sobre Γ_{1}, \\ \partial_{n} u + α u = g & sobre Γ_{2} . \end{cases} \end{matrix}

Formulación variacional

Observación: todos los cálculos que se hacen en esta sección son puramente formales, es decir, se supone siempre la regularidad suficiente para que todos ellos estén justificados.

$H^1(\Omega)$ $v\in L^2(\Omega)$ $L^2(\Omega)$ :

\begin{matrix} H^{1} (Ω) = {v | v \in L^{2} (Ω), \partial_{x} v, \partial_{y} v \in L^{2} (Ω)} \end{matrix}

$H^1(\Omega)$ $H^1(\Omega)$ $\Gamma_1$ :

\begin{matrix} V = {v | v \in H^{1} (Ω), v |_{Γ_{1}} = 0} . \end{matrix}

$\eqref{eq_P}$ se procede como sigue:

$\eqref{eq_P}$ $v \in V$ $\Omega$ :
$\begin{matrix} (1) & \int_{Ω} (- Δ u v + u v) d x d y = \int_{Ω} f v d x d y \forall v \in V . \end{matrix}$
Se aplica la fórmula de integración por partes:
$\begin{matrix} \int_{Ω} - Δ f g d x d y = \int_{Ω} \nabla f \nabla g d x d y - \int_{\partial Ω} (\partial_{n} f) g d σ . \end{matrix}$
$\eqref{eq_1}$
$\begin{matrix} \int_{Ω} \nabla u \nabla v d x d y + \int_{Ω} u v d x d y - \int_{\partial Ω} (\partial_{n} u) v d σ = \int_{Ω} f v d x d y \forall v \in V, \end{matrix}$
donde la integral de frontera se puede escribir
$\begin{matrix} - \int_{\partial Ω} (\partial_{n} u) v d σ = \int_{Γ_{2}} (\partial_{n} u) v d σ = \int_{Γ_{2}} (α u - g) v d σ . \end{matrix}$
$\eqref{eq_P}$ :
$\begin{matrix} (PV) & {\begin{cases} Hallar u \in {\tilde{u}}_{0} + V tal que \\ \int_{Ω} \nabla u \nabla v d x d y + \int_{Ω} u v d x d y + \int_{Γ_{2}} α u v d σ = \int_{Ω} f v d x d y + \int_{Γ_{2}} g v d σ \forall v \in V, \end{cases} \end{matrix}$
$\widetilde{u}_0$ $H^1(\Omega)$ $u_0$ $\Gamma_1$ $\widetilde{u}_0 \in H^1(\Omega)$ $\widetilde{u}_0 \big|_{\Gamma_1} = u_0$ .
$\eqref{eq_PV}$ es un problema de dimensión infinita: buscamos la solución en un espacio de Hilbert de dimensión infinita. Para abordar su resolución numérica necesitamos discretizarlo, esto es, sustituirlo por otro aproximado que tenga dimensión finita.

Discretización del problema variacional (PV) mediante Elementos Finitos P_1-Lagrange

$\Omega$ $\partial\Omega$ , está formada por una o varias curvas poligonales. Si no fuera así, el primer paso en la discretización sería aproximar el dominio por otro con dichas características.

$\bar\Omega$ $\cal T$ ${\cal T} = \{T_j\}_{j=1}^{R}$ verificando

\begin{matrix} \begin{array}{l} ⋃_{j = 1}^{R} T_{j} = \bar{Ω} \\ {\overset{˚}{T}}_{j} \neq \emptyset \forall j = 1, \dots, R \\ T_{i} \cap T_{j} = {\begin{cases} un lado com \overset{´}{u} n a ambos tri \overset{´}{a} ngulos o \\ un v \overset{´}{e} rtice com \overset{´}{u} n a ambos tri \overset{´}{a} ngulos o \\ \emptyset \end{cases} \forall i \neq j \end{array} \end{matrix}

Consideremos, por ejemplo, la de la Figura siguiente:

La bondad de la aproximación que vamos a describir a continuación dependerá, entre otras cosas, del tamaño de los triángulos que consideremos. Para medir éste se introduce el diámetro de la triangulación:

\begin{matrix} h = max_{T \in T} h_{T}, donde h_{T} = max_{x, y \in T} ∥ x - y ∥ \end{matrix}

${\cal T}_h$ $\bar\Omega$ $h$ .

Denotamos

\begin{matrix} W_{h} = {v_{h} | v_{h} \in C^{0} (\bar{Ω}), v_{h} |_{T} \in P_{1} (T) \forall T \in T_{h}} \end{matrix}

\begin{matrix} V_{h} = {v_{h} | v_{h} \in W_{h}, v_{h} |_{Γ_{1}} = 0} . \end{matrix}

$W_h$ $H^1(\Omega)$ $\text{dim\,} W_h = M$ $M$ $\mathcal{T}_h$ $V_h$ $V$ $\text{dim\,} V_h = M_0$ $M_0$ $\mathcal{T}_h$ $\Gamma_1$ .

$\{a_i\}_{i\in I}$ $\mathcal{T}_h$ $\{a_i\}_{i\in I_0}$ $I_0\subset I$ $\mathcal{T}_h$ $\Gamma_1$ $\text{card}(I) = M$ $\text{card}(I_0) = M_0$ .

$v_h\in W_h$ $\{v_h(a_i)\}_{i\in I}$ $v_h\in V_h$ $\{v_h(a_i)\}_{i\in I_0}$ .

$W_h$ $\{\varphi_i\}_{i\in I}$ definidas por

\begin{matrix} φ_{i} \in W_{h} \forall i = 1, \dots, M, φ_{i} (a_{j}) = δ_{i j} \forall i, j = 1, \dots, M . \end{matrix}

$\varphi_i$ $a_i$ $\mathcal{T}_h$ y solo es distinta de cero en los triángulos que contienen a dicho vértice.

$V_h$ $\{\varphi_i\}_{i\in I_0}$ $W_h$ . Luego,

\begin{matrix} v_{h} \in W_{h} ⟺ v_{h} (x) = \sum_{i \in I} v_{h} (a_{i}) φ_{i} (x), \end{matrix}

\begin{matrix} v_{h} \in V_{h} ⟺ v_{h} (x) = \sum_{i \in I_{0}} v_{h} (a_{i}) φ_{i} (x) . \end{matrix}

$\eqref{eq_PV}$ es la que sigue:

\begin{matrix} (PV_{h}) & {\begin{cases} Hallar u_{h} \in {\tilde{u}}_{0 h} + V_{h} tal que \\ \int_{Ω} \nabla u_{h} \nabla v_{h} d x d y + \int_{Ω} u_{h} v_{h} d x d y + \int_{Γ_{2}} α u_{h} v_{h} d σ = \int_{Ω} f v_{h} d x d y + \int_{Γ_{2}} g v_{h} d σ \\ \forall v_{h} \in V_{h}, \end{cases} \end{matrix}

$\widetilde{u}_{0h} \in W_h$ $\widetilde{u}_{0h}(a_i) = \widetilde{u}_{0}(a_i)$ $a_i\in \Gamma_1$ .

$\{\varphi_i\}_{i\in I_0}$ $V_h$ , el problema anterior se puede escribir

\begin{matrix} (2) & {\begin{cases} Hallar u_{h} \in {\tilde{u}}_{0 h} + V_{h} tal que \\ \int_{Ω} \nabla u_{h} \nabla φ_{i} d x d y + \int_{Ω} u_{h} φ_{i} d x d y + \int_{Γ_{2}} α u_{h} φ_{i} d σ = \int_{Ω} f φ_{i} d x d y + \int_{Γ_{2}} g φ_{i} d σ \forall i \in I_{0} . \end{cases} \end{matrix}

$u_h$ $W_h$ :

\begin{matrix} u_{h} \in {\tilde{u}}_{0 h} + V_{h} ⟺ u_{h} = \sum_{j \in I} u_{h} (a_{j}) φ_{j} = \sum_{j \in I_{0}} u_{h} (a_{j}) φ_{j} + \sum_{j \in I ∖ I_{0}} {\tilde{u}}_{0} (a_{j}) φ_{j}, \end{matrix}

$u_h$ $u_h(a_j)$ $j\in I_0$ .

$u_h$ $\eqref{eq_2}$ , el problema que tenemos que resolver se escribe:

\begin{matrix} (3) & {\begin{cases} Hallar u_{h} (a_{j})_{j \in I_{0}} tales que \\ \sum_{j \in I_{0}} u_{h} (a_{j}) (\int_{Ω} \nabla φ_{j} \nabla φ_{i} d x d y + \int_{Ω} φ_{j} φ_{i} d x d y + \int_{Γ_{2}} α φ_{j} φ_{i} d σ) \\ = \int_{Ω} f φ_{i} d x d y + \int_{Γ_{2}} g φ_{i} d σ \\ - \sum_{j \in I ∖ I_{0}} {\tilde{u}}_{0} (a_{j}) (\int_{Ω} \nabla φ_{j} \nabla φ_{i} d x d y + \int_{Ω} φ_{j} φ_{i} d x d y + \int_{Γ_{2}} α φ_{j} φ_{i} d σ) \\ \forall i \in I_{0} . \end{cases} \end{matrix}

Utilizando la siguiente notación

\begin{matrix} \begin{matrix} {\begin{cases} u_{j} = u_{h} (a_{j}), \\ a_{i j} = \int_{Ω} \nabla φ_{j} \nabla φ_{i} d x d y + \int_{Ω} φ_{j} φ_{i} d x d y + \int_{Γ_{2}} α φ_{j} φ_{i} d σ, \\ s_{i} = \int_{Ω} f φ_{i} d x d y + \int_{Γ_{2}} g φ_{i} d σ, \\ z_{i} = - \sum_{j \in I ∖ I_{0}} {\tilde{u}}_{0} (a_{j}) (\int_{Ω} \nabla φ_{j} \nabla φ_{i} d x d y + \int_{Ω} φ_{j} φ_{i} d x d y + \int_{Γ_{2}} α φ_{j} φ_{i} d σ), \end{cases} \end{matrix} \end{matrix}

$\eqref{eq_3}$ se puede reescribir de la forma, más compacta:

\begin{matrix} (4) & {\begin{cases} Hallar (u_{j})_{j \in I_{0}} tales que \\ \sum_{j \in I_{0}} a_{i j} u_{j} = s_{i} - z_{i} \forall i \in I_{0} \end{cases} \end{matrix}

$\eqref{eq_3}$ $M_0$ ecuaciones.

Bloqueo de las condiciones de contorno

$(a_{ij})_{i,j\in I_0}$ $(a_{ij})_{i,j\in I}$ .

$\eqref{eq_4}$ , sino el que vemos a continuación, obtenido penalizando el anterior.

$\eqref{eq_4}$ es equivalente a

\begin{matrix} (5) & {\begin{cases} Hallar (u_{j})_{j \in I} tales que \\ {\begin{cases} \sum_{j \in I_{0}} a_{i j} u_{j} = s_{i} & i \in I_{0} \\ u_{i} = {\tilde{u}}_{0} (a_{i}) & i \in I ∖ I_{0} . \end{cases} \end{cases} \end{matrix}

$u_h$ $\Gamma_1$ $u_0$ $M$ ecuaciones.

$\eqref{eq_5}$ es equivalente a

\begin{matrix} (6) & \sum_{j \in I} {\tilde{a}}_{i j} u_{j} = {\tilde{s}}_{i}, i \in I \end{matrix}

donde

\begin{matrix} (7) & \begin{matrix} {\tilde{a}}_{i j} = {\begin{cases} V & si i = j \in I ∖ I_{0}, \\ a_{i j} & si no, \end{cases} {\tilde{s}}_{i} = {\begin{cases} V {\tilde{u}}_{0} (a_{i}) & si i = j \in I ∖ I_{0}, \\ s_{i} & si no, \end{cases} \end{matrix} \end{matrix}

$V$ muy grande $\approx 10^{30}$ ).

$i\in I\setminus I_0$ $i$ $\eqref{eq_6}$ es

\begin{matrix} \sum_{j \in I} {\tilde{a}}_{i j} u_{j} = \sum_{j \in I, j \neq i} a_{i j} u_{j} + {\tilde{a}}_{i i} u_{i} = {\tilde{s}}_{i}, \end{matrix}

$u_i$ :

\begin{matrix} u_{i} = \frac{{\tilde{s}}_{i} - \sum_{j \in I, j \neq i} a_{i j} u_{j}}{{\tilde{a}}_{i i}} = \frac{V {\tilde{u}}_{0} (a_{i}) - \sum_{j \in I, j \neq i} a_{i j} u_{j}}{V} \approx \frac{V {\tilde{u}}_{0} (a_{i})}{V} \approx {\tilde{u}}_{0} (a_{i}) . \end{matrix}

$\eqref{eq_6}$ $i\in I_0$ $\eqref{eq_4}$ .

$\eqref{eq_6} + \eqref{eq_7}$ es el sistema que, realmente, se resuelve. Inicialmente, se construyen la matriz y el segundo miembro del sistema "sin penalizar"

\begin{matrix} \sum_{j \in I} a_{i j} u_{j} = s_{i}, i \in I \end{matrix}

$\eqref{eq_7}$ .

Cálculo de matrices y segundos miembros

$A$ $A=(a_{ij})$ $1\le i,j\le M$ $M$ $\mathcal{T}_h$ y

\begin{matrix} a_{i j} = \int_{Ω} \nabla φ_{i} \nabla φ_{j} d x d y + \int_{Ω} φ_{i} φ_{j} d x d y + \int_{Γ_{2}} α φ_{i} φ_{j} d σ, 1 \leq i, j \leq M . \end{matrix}

La condición de Dirichelet se introducirá más tarde, mediante la técnica del bloqueo ya explicada. Sea también

\begin{matrix} s_{i} = \int_{Ω} f φ_{i} d x d y + \int_{Γ_{2}} g φ_{i} d σ, 1 \leq i \leq M . \end{matrix}

$A$ es simétrica, definida positiva y hueca.

Nota: Se denomina matriz hueca aquélla cuyas componentes son mayoritariamente nulas. Las matrices de grandes dimensiones que tienen esta característica se suelen almacenar de formas especiales, no solo para evitar el almacenamiento de cantidades ingentes de ceros, sino también para evitar operaciones inútiles con ceros (multiplicaciones y sumas).

Una de las maneras de almacenar una matriz hueca es hacerlo como matriz sparse. Para ello se utilizan tres vectores: uno que contiene las componentes de la matriz que no son nulas y otros dos que contienen, respectivamente, los números de la fila y la columna en la que se encuentran.

Nota: Se llama matriz banda a una matriz que es hueca (es decir, cuyos elementos son mayoritariamente nulos) y cuyos elementos no nulos están concentrados en una banda de anchura pequeña (en relación con la dimensión de la matriz) en torno a la diagonal principal. Una manera de almacenar este tipo de matrices es el almacenamiento en orden perfil, consistente en guardar solamente las componentes de la matriz que están dentro de la mencionada banda.

$A$ se puede escribir como suma de tres matrices:

\begin{matrix} A_{r} = {(\int_{Ω} \nabla φ_{i} \nabla φ_{j} d x d y)}_{1 \leq i, j \leq M} A_{m} = {(\int_{Ω} φ_{i} φ_{j} d x d y)}_{1 \leq i, j \leq M} \end{matrix}

\begin{matrix} A_{Γ} = {(\int_{Γ_{2}} α φ_{i} φ_{j} d σ)}_{1 \leq i, j \leq M} \end{matrix}

Las dos primeras se llaman matriz de rigidez y matriz de masa, respectivamente.

Veamos cómo se realiza el cálculo de la matriz de rigidez, siendo análogo el cálculo para las otras dos. Denotamos

\begin{matrix} (8) & a_{i j}^{r} = \int_{Ω} \nabla φ_{i} \nabla φ_{j} d x d y = \sum_{T \in T_{h}} \int_{T} \nabla φ_{i} \nabla φ_{j} d x d y 1 \leq i, j \leq M . \end{matrix}

$\displaystyle\int_{T} \nabla \varphi_i \nabla \varphi_j \, dx$ $\varphi_i$ $\varphi_j$ $T$ $a_i$ $a_j$ $T$ $\eqref{eq_8}$ se reduce en realidad a unos pocos términos no nulos:

\begin{matrix} a_{i j}^{r} = \sum_{T ∋ a_{i}, a_{j}} \int_{T} \nabla φ_{i} \nabla φ_{j} d x d y . \end{matrix}

En la práctica, el cálculo de los elementos de la matriz no se hace elemento a elemento, ya que esto implicaría:

$a_{ij}$ $i$ $j$ ,
repetición de muchos cálculos.

$T$ $\displaystyle\int_{T} \nabla \varphi_i \nabla \varphi_j \, dx\,dy$ que no sean nulas) y acumularlos en los elementos correspondientes de la matriz. Obsérvese que, en principio, habría que calcular 9 valores:

\begin{matrix} \int_{T} \nabla φ_{i} \nabla φ_{j} d x d y i, j = 1, 2, 3. \end{matrix}

Pero, dado que

\begin{matrix} \int_{T} \nabla φ_{i} \nabla φ_{j} d x d y = \int_{T} \nabla φ_{j} \nabla φ_{i} d x d y, \end{matrix}

este número se reduce a 6. Esto queda descrito en el siguiente algoritmo.

Algoritmo 1 (cálculo de la matriz de rigidez)

Inicializar a cero los elementos de la matriz
$T_j$ $j = 1,2\dots, \text{nt}$ $i_1$ $i_2$ $i_3$ los números de sus tres vértices, calcular las 6 integrales
$\begin{matrix} (9) & \int_{T} | \nabla φ_{i_{1}} |^{2} d x d y, \int_{T} | \nabla φ_{i_{2}} |^{2} d x d y, \int_{T} | \nabla φ_{i_{3}} |^{2} d x d y, \end{matrix}$ $\begin{matrix} (10) & \int_{T} \nabla φ_{i_{1}} \nabla φ_{i_{2}} d x d y, \int_{T} \nabla φ_{i_{1}} \nabla φ_{i_{3}} d x d y, \int_{T} \nabla φ_{i_{2}} \nabla φ_{i_{3}} d x d y \end{matrix}$
$a^r_{ij}$ correspondiente, es decir, hacer
$\begin{matrix} a_{i_{1} i_{1}}^{r} ⟵ a_{i_{1} i_{1}}^{r} + \int_{T} | \nabla φ_{i_{1}} |^{2} d x d y (análogamente para las otras dos integrales en (9)) \end{matrix}$ $\begin{matrix} a_{i_{1} i_{2}}^{r} ⟵ a_{i_{1} i_{2}}^{r} + \int_{T} \nabla φ_{i_{1}} \nabla φ_{i_{2}} d x d y (análogamente para las otras dos integrales en (10)) \end{matrix}$

Detalle de los cálculos elementales

En construcción.

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla