θ(i,ξ,X) = prod([
(ξ - X[j])/(X[i+1] - X[j])
for j in setdiff(1:length(X),i+1)
] )θ (generic function with 1 method)Elementos Finitos de orden superior
Métodos Numéricos Avanzados. 2023-2024
David Gómez-Castro 
Más allá de simplexes y funciones lineales
¿Por qué simplexes y funciones lineales?
Hasta ahora, hemos trabajado con simplexes (segmentos, triángulos, …) y con funciones lineales en cada simplex.
La malla podría no de simplexes. Por ejemplo, en dimensión 2, podría ser de cuadrados.
Pero si \(K\) es un cuadrado de vértices \(\mathbf x_i\), no podemos encontrar una función lineal tal que \(\varphi(\mathbf x_j) = \delta_{1j}\).
Pero podemos encontrar \(\varphi\) polinomios de orden superior.
Elemento de referencia
Vamos a tomar un elemento “normalizado”, y a utilizarlo para construir las funciones de base. Digamos que la malla \(\mathcal T\) se compone de subconjunto \(K_i \in \mathbb T\), todos ellos transformaciones de uno de conjunto de referencia \(K_a\) (si hay varios \(K_a\) se habla de mallas mixtas).
La idea es construir un conjunto de polinomios \(P\) definidos \(K_a\), que se pueda “pegar” para construir una base \[ V_h = \{ v \in V : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \] donde \(\mathbf x_{K_i} : K_a \to K_i\).
En [Ciarlet, 2002] se hace estudia la noción rigurosa de “elemento”, y los elementos válidos y no válidos.
Polinomios de orden \(p > 1\) y conformes \(H^1\)
En esta sección buscamos discretizar \[ V_h = \{ v \in H^1(\Omega) : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \]
donde \(P\) son polinomios de grado \(p\).
Dado que en el interior de cada \(K_i\) tenemos polinomios, basta pegar de manera continua en las aristas.
Funciones nodales
Tomamos un elemento de referencia \(K_a = (-1,1)\). Vamos a construir polinomios de orden \(p\) que se anulen en ciertos puntos \[ X_0 = -1, \qquad X_1 = 1, \qquad X_{j} = -1 + \frac{2(j-1)}{p} \text{para} j = 2, \cdots, p \]
Así \[ \theta_i (\xi) = \frac{ \xi - X_1 }{X_i - X_1} \frac{ \xi - X_2 }{X_i - X_2} \cdots \frac{\xi - X_{i-1}} {X_i - X_{i-1}} \frac{\xi - X_{i+1}} {X_i - X_{i+1}} \cdots \frac{\xi - X_{p+1}} {X_i - X_{p+1}} \]
p = 2
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
using Plots, LaTeXStrings
plot(xlabel=L"\xi")
for i=0:p
plot!(ξ -> θ(i,ξ,X), label="i=$i", xlim=(-1,1))
end
scatter!(X, zeros(size(X)), label="")
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
plot(xlabel=L"\xi")
for i=0:p
plot!(ξ->θ(i,ξ,X), label="i=$i", xlim=(-1,1))
end
scatter!(X, zeros(size(X)), label="")
Funciones jerárquicas
La idea de las familias jerárquicas es encontrar bases \(\mathcal B_p \subset \mathcal B_{p+1}\), donde \(\mathcal B_p\) es una base de \(\mathbb R_p [X]\).
Para buscar bases de orden superior, se toman los polinomios de Legendre:
\[ L_n (x) = \frac{(-1)^n}{2^n n!} \frac{d^n}{dx^n} \left[ (1-x)^n(1+x)^n \right] \]
Tienen la propiedad
\[ \begin{aligned} L_0(x) &= 1, L_1(x) = x, \\ L_k(x) &= \frac{2k-1}{k} x L_{k-1} (x) - \frac{k-1}{k} L_{k-2}(x), \forall k \ge 2. \end{aligned} \]
y los polinomios de Lobatto
\[ \begin{aligned} l_0(x) &= \frac{1-x}{2}, \\ l_1(x) &= \frac{x+1}{2}, \\ l_k(x) &= \sqrt{\frac{2k-1}{2}} \int_{-1}^x L_{k-1} (\xi) d \xi. \end{aligned} \]
Las funciones \(l_0, \cdots, l_p\) forman una base de \(\mathbb R_p[X]\).
Además \(l_k(-1) = l_k(1) = 0\) para \(k \ge 2\).
En este contexto \(\theta_k = l_k\) para \(k = 0, \cdots, p\).
Construcción de las funciones de base
Fijamos un intervalo \((a,b)\) y una malla \(K_i = (x_{i}, x_{i+1})\) con \(i=1, \cdots I\) donde es construir la aplicación afín \(\mathbf x_i : (-1,1) \to (x_{i}, x_{i+1})\).
Tomemos \(p\) fijo y una elección de \(\theta_i\) una de las elecciones anteriores (nodal o jerárquica). Observamos que \[ \begin{aligned} \theta_0 (-1) = 1 &, \qquad \theta_0(1) = 0, \\ \theta_1 (-1) = 0 &, \qquad \theta_1(1) = 1, \\ \theta_k(-1) = 0 &, \qquad \theta_k(1) = 0. \end{aligned} \]
Para los problemas elípticos que hemos visto, queremos que la base para Galerkin sea de funciones continuas tenemos que pegar con continuidad.
Esto nos permite construir funciones base de tipo vértice (con \(v_{i,0} (x_j) = \delta_{ij}\)) \[ v_{i,0}(x) = \begin{cases} (\theta_1 \circ \mathbf x_{K_i}^{-1}) (x) , & x \in K_i \\ (\theta_0 \circ \mathbf x_{K_{i+1}}^{-1}) (x) , & x \in K_{i+1} \\ 0, & x \notin K_{i} \cup K_{i+1}. \end{cases} \]
Ejemplo. Tomando como \(\theta_i\) las funciones nodales
mesh = [0.1, 0.2, 0.4]
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
function v(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(1,ξ,X)
elseif x >= mesh[2] && x < mesh[3]
ξ = -1 + 2*(x-mesh[2])/(mesh[3] - mesh[2])
return θ(0,ξ,X)
else
return 0.0
end
end
plot(x -> v(x), label=L"v", xlim=(0,0.5), xlabel=L"x")
scatter!(mesh, zeros(size(mesh)), label="")
También podemos construir funciones base de tipo burbuja, para \(k \ge 2\)
\[ v_{i,k-1}(x) = \begin{cases} \theta_{k} \circ \mathbf x_{K_i}^{-1} & x \in K_i \\ 0 & x \notin K_i \end{cases} \]
Estas no existen para \(p = 1\).
Ejemplo. Tomando como \(\theta_i\) las funciones nodales
mesh = [0.1, 0.2, 0.4]
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
function w1(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(2,ξ,X)
else
return 0.0
end
end
function w2(x::Real)
if x >= mesh[1] && x < mesh[2]
ξ = -1 + 2*(x-mesh[1])/(mesh[2] - mesh[1])
return θ(3,ξ,X)
else
return 0.0
end
end
plot(x -> w1(x), label=L"w_1", xlim=(0,0.5))
plot!(x -> w2(x), label=L"w_2", xlim=(0,0.5))
scatter!(mesh, zeros(size(mesh)), label="")
xlabel!(L"x")
Espacio de funciones base
Con estas elecciones tomamos como espacio de funciones base para la aproximación de Galerkin \[ V_h = \mathrm{span} \{ v_{i,k} : i = 1, \cdots, I \text{ y } k = 0, \cdots, p \}. \]
Este mismo planteamiento se puede hacer en dimensión superior, cuando \(K\) no es un triángulo/tetrahedro: [Šolin, Segeth & Doležel, 2004]
Usando SymPy.jl
using SymPy, Plots, LaTeXStrings
Δx = 1.0
xx = -1.0:Δx:1.0
ξ = symbols("ξ", real=true)
x = symbols("x", real=true)
p = 3
X = [-1; 1; [-1 + 2*j/p for j=1:p-1]]
θ(i,ξ,X) = prod([
(ξ - X[j])/(X[i+1] - X[j])
for j in setdiff(1:length(X),i+1)
] )
v(i,x) = sympy.Piecewise(
( θ(1,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
( θ(0,-1 + 2*(x-xx[i])/Δx, X), And(Ge(x,xx[i]),Le(x, xx[i+1]))) ,
(0, True))
w1(i,x) = sympy.Piecewise(
( θ(2,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
(0, True))
w2(i,x) = sympy.Piecewise(
( θ(3,-1 + 2*(x-xx[i-1])/Δx, X), And(Ge(x,xx[i-1]),Le(x, xx[i]))) ,
(0, True))
φ = [ [ v(i,x) for i in 2:length(xx)-1 ];
[ w1(i,x) for i in 2:length(xx) ] ;
[ w2(i,x) for i in 2:length(xx) ] ;
];
plot(legend_position=:outerright,xlim=(-1.01,1.01))
scatter!(xx,zeros(length(xx)),label="")
for i in 1:length(φ)
plot!(φ[i],label=latexstring("\\varphi_{$i}"))
end
plot!()
xlabel!(L"x")
Dφ = diff.(φ,x)
NN = length(φ)
A = zeros(NN,NN)
for i=1:NN, j=1:NN
A[i,j] = integrate(Dφ[j]*Dφ[i],(x,-1.0,1.0))
end
A
u_analytic = 1 - x^4
f_analytic = 12*x^2
b = [integrate(f_analytic*φ[i],(x,-1.0,1.0)) for i=1:NN]
uh_coord = A \ b;
uh(x) = sum([uh_coord[i]*φ[i](x) for i=1:length(φ)])
plot(x->uh(x), label="Sol numérica", thickness_scaling=1.5,linewidth=2,legend=:outerright,xlim=(-1.0001,1.0001))
plot!(u_analytic, label="Sol exacta",linewidth=2)
scatter!(xx,zeros(length(xx)),label="Malla")
xlabel!(L"x")
ylabel!(L"u")
Funciones de base con mayor regularidad
Problemas de orden superior
Las construcciones anteriores son solo continuas (no \(C^1\)) en los puntos de la malla. Esto es suficiente para problemas de segundo orden como los que hemos visto: \(-\frac{\partial^2 u}{\partial x^2} = f\).
Es muy habitual, por ejemplo en elasticidad, encontrar problemas de cuarto ordern \(\frac{\partial^4 u}{\partial x^4} = f\). La formulación débil es, entonces (si \(I = (-1,1)\)) \[ \int_{-1}^1 \frac{\partial^2 u}{\partial x^2} \frac{\partial^2 \varphi}{\partial x^2} = \int_{-1}^1 f \varphi + \text{ términos en $-1$ y $1$ (condiciones de contorno)} \]
En este caso buscaremos discretizar \[ V_h = \{ v \in H^2(\Omega) : v|_{K_i} \circ \mathbf x_{K_i} \in P \} \]
Debemos “pegar” de manera \(C^1\) en las aristas de la malla.
Elemento de Hermite en \([0,1]\)
Para “pegar” de forma \(C^1\) debemos tener cuidado con los valores de \(\varphi'\) en los puntos de la malla.
Si trabajamos en el intervalo de referencia \(K = [0,1]\) esto significa tener cuidado con \(\varphi'(0), \varphi'(1)\).
Podemos construir primero dos polinomios de base cuyas derivadas se anulen en los extremos ellas o sus derivadas \[ \theta_0(0) = 1, \qquad \theta_0(1) = 0, \qquad \theta_0'(0) = 0, \qquad \theta_0'(1) = 0 \]
\[ \theta_1(0) = 0, \qquad \theta_1(1) = 1, \qquad \theta_1'(0) = 0, \qquad \theta_1'(1) = 0 \]
\[ \theta_2(0) = 0, \qquad \theta_2(1) = 0, \qquad \theta_2'(0) =1, \qquad \theta_2'(1) = 0 \]
\[ \theta_3(0) = 0, \qquad \theta_3(1) = 0, \qquad \theta_3'(0) =0, \qquad \theta_3'(1) = 1 \]
Esto requiere, al menos, grado 3, y en tal caso
θHermite = [
ξ -> 2ξ^3 - 3ξ^2 + 1, # θHermite[i] = θ_{i-1}
ξ -> -2ξ^3 + 3ξ^2,
ξ -> ξ^3 - 2ξ^2 + ξ,
ξ -> ξ^3 - ξ^2
]
plot(xlim=(0,1),xlabel=L"\xi")
for i=1:4
plot!(ξ->θHermite[i](ξ),label=latexstring("\\theta_{$(i-1)}"))
end
plot!()
Ahora pegamos haciendo \(v_{i}(x_j) = \delta_{ij}, v'_i (x_j) = 0\) y \(w_i (x_j) = 0, w_i'(x_j) = \delta_{ij}\) \[ v_{i} (x) = \begin{cases} \theta_1 \left( \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \right) , & x \in [x_{i-1},x_{i}) \\ \theta_0 (\frac{x-x_{i}}{x_{i+1}-x_{i}}) , & x \in x \in [x_{i},x_{i+1})\\ 0, & \text{otherwise}. \end{cases} \qquad w_{i} (x) = \begin{cases} (x_{i}-x_{i-1}) \theta_4 \left( \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \right) , & x \in [x_{i-1},x_{i}) \\ (x_{i+1}-x_{i}) \theta_3 (\frac{x-x_{i}}{x_{i+1}-x_{i}}) , & x \in x \in [x_{i},x_{i+1})\\ 0, & \text{otherwise}. \end{cases} \]
mesh = [-Inf, 0.0, 0.1, 0.2, 0.4, 0.5, +Inf]
function v(i::Integer,x::Real)
if x >= mesh[i-1] && x < mesh[i]
ξ = (x-mesh[i-1])/(mesh[i] - mesh[i-1])
return θHermite[2](ξ)
elseif x >= mesh[i] && x < mesh[i+1]
ξ = (x-mesh[i])/(mesh[i+1] - mesh[i])
return θHermite[1](ξ)
else
return 0.0
end
end
scatter(mesh[2:end-1], zeros(size(mesh[2:end-1])),
label="",
leg=:outerright)
for i=2:length(mesh)-1
plot!(x -> v(i,x), label=latexstring("v_{$(i-1)}"))
end
xlabel!(L"x")
mesh = [-Inf, 0.0, 0.1, 0.2, 0.4, 0.5, +Inf]
function v(i::Integer,x::Real)
if x >= mesh[i-1] && x < mesh[i]
dx = mesh[i] - mesh[i-1]
ξ = (x-mesh[i-1])/dx
return dx*θHermite[4](ξ)
elseif x >= mesh[i] && x < mesh[i+1]
dx = mesh[i+1] - mesh[i]
ξ = (x-mesh[i])/dx
return dx*θHermite[3](ξ)
else
return 0.0
end
end
scatter(mesh[2:end-1], zeros(size(mesh[2:end-1])),
label="",
leg=:outerright)
for i=2:length(mesh)-1
plot!(x -> v(i,x), label=latexstring("w_{$(i-1)}"))
end
xlabel!(L"x")
Bogner-Fox-Schmit Rectangular Element in 2D
Tomando \(K = [0,1] \times [0,1]\) podemos construir \[ \theta_{ij}(x,y) = \theta_i(x) \theta_j(y), \qquad i,j \in \{0,\cdots,3\} \]
x = y = 0:1e-1:1.0
for (i,j) in [(1,2),(4,1),(3,2),(3,3)]
f(x,y) = θHermite[i](x)*θHermite[j].(y)
z = f.(x',y) ; # note the ' which permutes the dims of x
display(surface(x, y, z,
title=latexstring("i=$(i-1), j=$(j-1)"),
xlabel=L"\xi_1",
ylabel=L"\xi_2"
))
end
También se pueden construir elementos de tipo \(C^1\) en el triángulo de referencia, e.g., triángulo de Argyris.
Para pegar estos elementos en una malla hay que tener cuidado con rescalar las orientación de los bordes (ver [Domínguez & Sayas, 2008], [Valdman, 2020]).
Más información [Ciarlet & Lions, 1990, Capítulo VII]
References
Ciarlet, P.G. (2002) Re-publishment of the 1978 edition. The Finite Element Method for Elliptic Problems. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898719208.
Ciarlet, P.G. & Lions, J.-L. (1990) Handbook of numerical analysis. Amsterdam ; New York : New York, North-Holland ; Distributors for the United States and Canada, Elsevier Science Pub. Co.
Domínguez, V. & Sayas, F.-J. (2008) Algorithm 884: A Simple Matlab Implementation of the Argyris Element. ACM Transactions on Mathematical Software. 35 (2), 1–11. doi:10.1145/1377612.1377620.
Šolin, P., Segeth, K. & Doležel, I. (2004) Higher-order finite element methods. Studies in advanced mathematics. Boca Raton, FL, Chapman & Hall/CRC.
Valdman, J. (2020) MATLAB Implementation of C1 Finite Elements: Bogner-Fox-Schmit Rectangle. In: R. Wyrzykowski, E. Deelman, J. Dongarra, & K. Karczewski (eds.). Parallel Processing and Applied Mathematics. Cham, Springer International Publishing. pp. 256–266. doi:10.1007/978-3-030-43222-5_22.