Elementos Finitos de orden superior

Métodos Numéricos Avanzados. 2023-2024

David Gómez-Castro

Universidad Complutense de Madrid

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}} \]

θ(i,ξ,X) = prod([ 
- X[j])/(X[i+1] - X[j]) 
        for j in setdiff(1:length(X),i+1) 
        ] )
θ (generic function with 1 method)
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")

Usando SymPy.jl

= 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