¿Por qué simplexes y funciones lineales?

Hasta ahora, hemos trabajado con simplexes (segmentos, triángulos, …) y con funciones lineales en cada simplex.

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).

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)

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})\).

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

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\)

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

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")

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\).

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.

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

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]