\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2019 26(2) : 215–251 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v26i2. 35968

una implementación computacional del

método vem mixto para el problema de brinkman en 2d

a computational implementation of the mixed-vem method for the brinkman problem in 2d

Helen Guillén-Oviedo^{∗}  Filánder A. Sequeira^{†}

Received: 19/Feb/2019; Revised: 9/May/2019; Accepted: 23/May/2019

∗Universidad Nacional, Escuela de Matemática, Heredia, Costa Rica. E-Mail: hellen.guillen.oviedo@una.ac.cr

†Misma dirección que/Same address as: H. Guillén-Oviedo. E-Mail: filander.sequeira.chavarria@una.ac.cr

216 h. guillén — f.a. sequeira

Resumen

En este artículo se describen algunos aspectos específicos sobre una implementación computacional para la formulación mixta de elementos virtuales (mixed-VEM, por sus siglas en inglés) del problema lineal de Brinkman en dos dimensiones, con condiciones de frontera de Dirichlet no homogéneas. La formulación empleada fue originalmente propuesta y analizada en [18]. La implementación planteada aquí considera cualquier grado polinomial k ≥ 0 de manera natural al construir diversas matrices locales de bajo tamaño. Además, se propone un algoritmo para el ensamblaje del sistema lineal global asociado, que garantiza la continuidad de la componente normal en la solución discreta.

Palabras clave: modelo de Brinkman; método mixto de elementos virtuales; aproximaciones de alto orden; implementación computacional.

Abstract

In this paper we describe some specific aspects on the computational implementation of the a mixed virtual element method (mixed-VEM) for the two-dimensional linear Brinkman model with non-homogeneous Dirichlet boundary conditions. The formulation used below was originally proposed and analysed in [18]. The implementation presented here consider any polynomial degree k ≥ 0 in a natural way by building several local matrices of small size. In addition, an algorithm is proposed for the assembly of the associated global linear system, which guarantees the continuity of the normal component in the discrete solution.

Keywords: Brinkman model; mixed virtual element method; high-order approximations; computational implementation.

Mathematics Subject Classification: 65N30, 65K05.

1 Introducción

El problema de Brinkman describe un flujo que presenta algunos fluidos viscosos en un medio altamente poroso, y puede verse como una combinación dependiente de los parámetros de las ecuaciones de Darcy y Stokes. En particular, ha sido un modelo de estudio muy relevante para la comunidad de Análisis Numérico en las últimas décadas, dado que el mismo además puede asociarse a la discretización temporal del problema de Stokes evolutivo. Por ejemplo, algunos trabajos relacionados a este problema se presentan en [1, 2, 3, 23, 24, 27, 28, 33], así como las referencias en ellos.

En particular, una formulación dual-mixta con su respectivo esquema H(div)-conforme, combinado con el método de penalización interior simétrico, es presentada en [29]. Por otra parte, una formulación mixta en la cual el flujo de vorticidad se considera como una incógnita adicional fue introducida en [33]. Este enfoque difiere de otras formulaciones basadas en vorticidad, ya que utiliza la secuencia de Rham para derivar una ecuación para esta incógnita, en lugar de obtenerla al tomar el rotacional de la ecuación de momento.

Por otro lado, el método de elementos virtuales (VEM, por sus siglas en inglés), fue introducido en [6] para el problema modelo de Poisson, y es uno de los esquemas de discretización de alto orden que surgió como consecuencia natural de los nuevos desarrollos e interpretaciones del método mimético de diferencias finitas, ver por ejemplo, [13]. Las principales ventajas de los enfoques VEM incluyen una extensión de la técnica clásica de elementos finitos a las mallas poligonales y polihédricas generales, además de una generalización de los métodos miméticos a grados arbitrarios de precisión y propiedades de continuidad arbitrarias. Asimismo, tal y como se señala en [8], otro beneficio de los métodos VEM, consiste en la sólida base matemática para su análisis, en comparación con los métodos de volúmenes finitos. Luego, en el contexto de técnicas de elementos virtuales puramente mixtas, estas se desarrollaron inicialmente en [14], y más recientemente se extendió en [9, 10, 17]. En particular, [17] proporciona el primer análisis de un método de elementos virtuales para una formulación variacional mixta del problema de Stokes en el que el pseudo-esfuerzo y la velocidad son las únicas incógnitas, mientras que la presión se calcula mediante un posprocesamiento, ver también [16] para más detalles. Más aún, el análisis expuesto en [17] permite estudiar problemas de la misma naturaleza que Stokes, como ocurre con el problema de Brinkman, donde en [18] se introdujo y se analizó un método mixto de elementos virtuales. Para otras contribuciones sobre los métodos VEM, se considera, por ejemplo, a [4, 7, 12, 15, 19, 20, 22, 25, 26].

A diferencia de las referencias mencionadas previamente, el propósito del presente artículo consiste en describir una implementación computacional del método de elementos virtuales para aproximaciones de alto orden y sin estar ligada a un lenguaje de programación. Más precisamente, se detalla una forma de implementar el método mixto de elementos virtuales para el problema de Brinkman en dos dimensiones, cuyo esquema de Galerkin fue introducido y analizado en [18]. En particular, la codificación descrita aquí fue la empleada en [18], donde es posible verificar la calidad de esta por medio de los resultados numéricos proporcionados. La razón principal en esta contribución es que se cuenta con poca literatura que explique cómo programar los subespacios de elementos virtuales (en particular los mixtos) de manera precisa.

Algunos trabajos, como por ejemplo, [6, 8, 9, 12, 14, 16, 22], explican aspectos específicos sobre la implementación computacional, pero en general estos no describen las estructuras empleadas. Al menos por el conocimiento de los autores, solo en [5] se presenta un paradigma detallado de programación como el presentado más adelante, pero para una formulación primal del problema de elasticidad. Por ende, este artículo es el primero donde se detalla una implementación para elementos virtuales con subespacios ℍ(div)-conformes para el problema de Brinkman.

El artículo está organizado como sigue: en la Sección 2 se introduce el problema de valores de frontera de interés, y se recuerda su fomulación mixta basada en pseudo-esfuerzo. Luego, en la Sección 3 se resume el esquema de elementos virtuales propuesto en [18] para el problema modelo. Seguidamente, la contribución principal de este artículo se presenta en las secciones 4 y 5. En efecto, en la Sección 4 se detalla la forma de ensamblar el sistema lineal asociado al problema discreto de elementos virtuales mixto, mientras que en la Sección 5 se describen los aspectos computaciones requeridos para la construcción de las matrices locales que contribuyen a dicho ensamblaje. Es importante aclarar que todo esto se realiza para un grado polinomial k ≥ 0 arbitrario, lo que significa que la implementación obtenida es de alto orden. Finalmente, se presentan algunas conclusiones en la Sección 6.

Notaciones

Para finalizar esta sección se mencionan algunas notaciones a utilizar a lo largo de este escrito. En efecto, dado Ω ⊂ ℝ^{2} un dominio acotado con frontera Γ, se define n como el vector normal unitario exterior a Γ. Luego, se adopta la terminología estándar para los espacios de Lebesgue L^{p}(Ω), p > 1, y los espacios de Sobolev H^{s}(Ω), s ∈ ℝ, con norma ∥ · ∥_{s,Ω} y seminorma | · |_{s,Ω}. En particular, H^{1/2}(Γ) es el espacio de trazas de funciones de H^{1}(Ω), y H^{−1/2}(Γ) denota su respectivo dual. Además, por M y 𝕄 se denotan las contrapartes vectorial y tensorial del espacio funcional escalar arbitrario M. Así, denotando div (resp. rot) el operador usual de la divergencia div (resp. el operador rotacional rot) actuando sobre las filas de un tensor dado, se recuerdan los espacios:

H(div; Ω) := {τ ∈ L^{2}(Ω) : div(τ) ∈ L^{2}(Ω)} ,
ℍ(div; Ω) := {τ ∈ 𝕃^{2}(Ω) : div(τ) ∈ L^{2}(Ω)} ,
H(rot; Ω) := {τ ∈ L^{2}(Ω) : rot(τ) ∈ L^{2}(Ω)} ,

y

ℍ(rot; Ω) := {τ ∈ 𝕃^{2}(Ω) : rot(τ) ∈ L^{2}(Ω)} ,

equipados con las normas usuales

∥τ∥_{div;Ω}^{2} := ∥τ∥_{0,Ω}^{2} + ∥div(τ)∥_{0,Ω}^{2}, ∀ τ ∈ H(div; Ω) ,
∥τ ∥_{div;Ω}^{2} := ∥τ ∥_{0,Ω}^{2} + ∥div(τ)∥_{0,Ω}^{2}, ∀ τ ∈ ℍ(div; Ω) ,
∥τ∥_{rot;Ω}^{2} := ∥τ∥_{0,Ω}^{2} + ∥rot(τ)∥_{0,Ω}^{2}, ∀ τ ∈ H(rot; Ω) ,

y

∥τ ∥_{rot;Ω}^{2} := ∥τ ∥_{0,Ω}^{2} + ∥rot(τ)∥_{0,Ω}^{2}, ∀ τ ∈ ℍ(rot; Ω) ,

son espacios de Hilbert. Finalmente, dados τ := (τ_{ij}), ζ := (ζ_{ij}) ∈ ℝ^{2×2}, se escribe:

τ^{t} := (τ_{ji}) ,  tr(τ) := ∑_{i=1}^{2} τ_{ii} ,  τ^{d} := τ − \frac{1}{2} tr(τ) I ,

y

τ : ζ := ∑_{i,j=1}^{2} τ_{ij}ζ_{ij} ,

donde I es la matriz identidad de ℝ^{2×2}.

2 El problema de Brinkman y su formulación mixta

Sea Ω un dominio poligonal acotado y simplemente conexo de ℝ^{2} con frontera Γ. El objetivo es determinar la velocidad u, el pseudo-esfuerzo σ, y la presión p de un flujo constante de Brinkman que ocupa la región Ω.

En otras palabras, dado f ∈ L^{2}(Ω) una fuente y g ∈ H^{1/2}(Γ) el dato de Dirichlet, se busca un campo tensorial σ, un campo vectorial u y un campo escalar p tales que:

σ = μ ∇u − p I en Ω , α u − div(σ) = f en Ω ,
div(u) = 0 en Ω ,  u = g sobre Γ ,  ∫_{Ω} p = 0 ,

(1)

donde μ es la viscosidad dinámica, y α > 0 es una aproximación constante de la viscosidad dividida por la permeabilidad. Además, tal y como lo requiere la condición de incompresibilidad, se supone a partir de ahora que el dato g satisface la condición de compatibilidad ∫_{Γ} g · n = 0.

Por otro lado, la condición de incompresibilidad implica además que (1) puede ser reescrito de la forma:

\frac{1}{μ} σ^{d} = ∇u en Ω , α u − div(σ) = f en Ω ,
u = g sobre Γ , ∫_{Ω} tr(σ) = 0 ,

(2)

donde la presión p puede ser recuperada a través del posprocesamiento dado por:

p = −\frac{1}{2} tr(σ) en Ω .  (3)

De esta forma, procediendo tal y como en [23], se tiene que la velocidad es reemplazada por la segunda ecuación de (2), tal y como sigue:

u = \frac{1}{α} \{f + div(σ)\} en Ω ,  (4)

la cual da lugar a la siguiente formulación variacional dual-mixta de (2): hallar σ ∈ ℍ tal que:

a(σ, τ) = F (τ)  ∀ τ ∈ ℍ ,  (5)

donde:

ℍ := ℍ_{0}(div; Ω) := {τ ∈ ℍ(div; Ω) : ∫_{Ω} tr(τ) = 0} ,

a : ℍ × ℍ → ℝ es la forma bilineal:

a(ζ, τ) := \frac{1}{μ} ∫_{Ω} ζ^{d} : τ^{d} + \frac{1}{α} ∫_{Ω} div(ζ) · div(τ)  ∀ ζ, τ ∈ ℍ ,

y F : ℍ → ℝ es el funcional lineal:

F (τ) := −\frac{1}{α} ∫_{Ω} f · div(τ) + ⟨τ n, g⟩_{Γ}  ∀ τ ∈ ℍ .  (6)

La solubilidad única del esquema (5) es establecida en [18, Thm 2.1], y recordada en el siguiente teorema.

Teorema 2.1. Suponga que f ∈ L^{2}(Ω) y g ∈ H^{1/2}(Γ). Entonces, existe una única solución σ ∈ ℍ de (5). Además, existe C > 0 tal que:

∥σ∥_{div;Ω} ≤ C {∥f∥_{0,Ω} + ∥g∥_{1/2,Γ}}.

3 El esquema mixto de elementos virtuales

Para iniciar, considere {T_{h}}_{h>0} una familia de descomposiciones de Ω en elementos poligonales, donde h denota el mayor de sus diámetros. Por el bien de la teoría, se asume que T_{h} satisface las condiciones descritas al inicio de [18, Sección 3.1]. Ahora, de acuerdo con [18, Sección 5], dado un entero k ≥ 0, se considera el subespacio de elementos virtuales ℍ_{h} de ℍ := ℍ_{0}(div; Ω) dado por:

ℍ_{h} := {τ ∈ ℍ : τ |_{K} ∈ ℍ_{k}^{K} ∀ K ∈ T_{h}} ,  (7)

donde, para cada K ∈ T_{h}, se introduce el espacio local de elementos virtuales de orden k, ver por ejemplo, [9, 10]:

H_{k}^{K} := {τ := (τ_{1}, τ_{2})^{t} ∈ H(div; K) ∩ H(rot; K) : τ · n|_{e} ∈ P_{k}(e)
∀ arista e ∈ ∂K , div(τ) ∈ P_{k}(K) , y rot(τ) ∈ P_{k−1}(K)} ,

(8)

donde rot(τ) := \frac{∂τ_{2}}{∂x_{1}} − \frac{∂τ_{1}}{∂x_{2}} y P_{ℓ}(K) corresponde al espacio de polinomios definidos sobre K de grado total a lo más ℓ, con P_{−1}(K) := {0}. Más aún, se define para cada K ∈ T_{h} la versión tensorial de H_{k}^{K} dada por:

ℍ_{k}^{K} := {τ ∈ ℍ(div; K)∩ℍ(rot; K) : (τ_{i1}, τ_{i2})^{t} ∈ H_{k}^{K} ∀ i ∈ {1, 2}} , (9)

la cual se utiliza para definir ℍ_{h} en (7).

Seguidamente, se define el esquema de Galerkin asociado con (5) de la forma: hallar σ_{h} ∈ ℍ_{h} tal que:

a_{h}(σ_{h}, τ_{h}) = F (τ_{h})  ∀ τ_{h} ∈ ℍ_{h} ,  (10)

donde a_{h} : ℍ_{h} × ℍ_{h} → ℝ es la forma bilineal discreta global dada por:

a_{h}(ζ, τ) := ∑_{K∈T_{h}} a_{h}^{K}(ζ, τ)  ∀ ζ, τ ∈ ℍ_{h} ,  (11)

mientras que, para cada K ∈ T_{h}, a_{h}^{K} : ℍ_{h}^{K} × ℍ_{h}^{K} → ℝ es la forma bilineal discreta local:

a_{h}^{K}(ζ, τ) := \frac{1}{μ} ∫_{K} (Π̂_{k}^{K}(ζ))^{d} : (Π̂_{k}^{K}(τ))^{d} + \frac{1}{α} ∫_{K} div(ζ) · div(τ)
+ S^{K} (ζ − Π̂_{k}^{K}(ζ), τ − Π̂_{k}^{K}(τ)) ∀ ζ, τ ∈ ℍ_{h}^{K} .

(12)

La forma bilineal S^{K} está definida a través de los grados de libertad ℍ_{k}^{K}, mientras que Π̂_{k}^{K} es un operador de proyección adecuado sobre ℍ_{k}^{K}. Ambos operadores se describen en las siguientes secciones, respectivamente. Previo a esto, se comenta que el esquema discreto (10) posee solución única, lo cual fue probado en [18, Thm 5.1], mientras que en [18, Thm 5.3] se dedujeron las respectivas estimaciones de error a-priori.

3.1 Grados de libertad y la forma bilineal S^{K}

Sea e una arista de T_{h} con punto medio x_{e} y longitud h_{e}. Entonces, se considera el siguiente conjunto de k + 1 monomios normalizados sobre e:

B_{k}(e) := \{(\frac{x − x_{e}}{h_{e}} )^{j}\}_{0 ≤ j ≤ k} ,  (13)

el cual sin duda constituye una base de P_{k}(e). Similarmente, dado un entero ℓ ≥ 0 y un elemento K ∈ T_{h} con baricentro x_{K}, se define el siguiente conjunto de \frac{1}{2}(ℓ + 1)(ℓ + 2) monomios normalizados sobre K:

B_{ℓ}(K) := \{(\frac{x − x_{K}}{h_{K}} )^{α}\}_{0 ≤ |α| ≤ ℓ} ,  (14)

el cual es una base de P_{ℓ}(K). Es importante aclarar que en (14) se emplea la notación de multi-índices, donde, dado x := (x_{1}, x_{2})^{t} ∈ ℝ^{2} y α := (α_{1}, α_{2})^{t}, con enteros no negativos α_{1} y α_{2}, se tiene que x^{α} := x_{1}^{α_{1}}x_{2}^{α_{2}} y |α| := α_{1} + α_{2}. Seguidamente, se recuerdan de [9, 10] los siguientes grados de libertad locales para un τ ∈ H_{k}^{K}:

m_{q,e}^{n}(τ) := ∫_{e} τ · n q ∀ q ∈ B_{k}(e) , ∀ arista e ∈ ∂K ,
m_{q,K}^{div}(τ) := ∫_{K} τ · ∇q ∀ q ∈ B_{k}(K) ∖ {1} ,
m_{q,K}^{rot}(τ) := ∫_{K} τ · q ∀ q ∈ G_{k}^{⊥}(K) ,

(15)

donde G_{k}^{⊥}(K) es una base de (∇P_{k+1}(K))^{⊥} ∩ P_{k}(K), la cual corresponde al L^{2}(K)-ortogonal de ∇P_{k+1}(K) en P_{k}(K).

Así, gracias a las cardinalidades de B_{k}(e) y B_{k}(K), y de acuerdo con las dimensiones de P_{k}(K) y ∇P_{k+1}(K) (ver además [16, Lemma 4.3.1] por detalles), se obtiene que la cardinalidad de G_{k}^{⊥}(K) es \frac{1}{2}k(k + 1), y con ello, la cantidad de grados de libertad locales definidos en (15) viene dada por:

n_{k}^{K} := (k + 1) d_{K} + \{\frac{(k + 1)(k + 2)}{2} − 1\} + \frac{k(k + 1)}{2}
= (k + 1)(d_{K} + k + 1) − 1 ,  (16)

donde d_{K} corresponde al número de aristas en K. Más aún, en [9, Sección 3.4] se comprobó que, para cada K ∈ T_{h}, estos n_{k}^{K} grados de libertad locales son unisolventes en H_{k}^{K}.

De acuerdo a lo anterior, es natural considerar para (9), los siguientes 2n_{k}^{K} grados de libertad locales:

m_{q,e}^{n}(τ) := ∫_{e} τ n · q  ∀ q ∈ B_{k}(e) , ∀ arista e ∈ ∂K ,
m_{q,K}^{div}(τ) := ∫_{K} τ : ∇q  ∀ q ∈ B_{k}(K) ∖ {(1, 0)^{t}, (0, 1)^{t}} , (17)
m_{ρ,K}^{rot}(τ) := ∫_{K} τ : ρ  ∀ ρ ∈ G_{k}^{⊥}(K) ,

donde

B_{k}(e) := {(q, 0)^{t} : q ∈ B_{k}(e)} ∪ {(0, q)^{t} : q ∈ B_{k}(e)},
B_{k}(K) := {(q, 0)^{t} : q ∈ B_{k}(K)} ∪ {(0, q)^{t} : q ∈ B_{k}(K)} ,

y

G_{k}^{⊥}(K) := \{( \begin{matrix}q \\ 0 \end{matrix} ) : q ∈ G_{k}^{⊥}(K)\} ∪ \{( \begin{matrix}0 \\ q \end{matrix} ) : q ∈ G_{k}^{⊥}(K)\} .

En virtud de los grados de libertad definidos en (15) y (17), es posible definir una forma bilineal S^{K} : ℍ_{h}^{K} × ℍ_{h}^{K} → ℝ basada en estos. En efecto, sea K ∈ T_{h} y considere la unión de todos los K-momentos (ver (15)) de un determinado τ ∈ H^{1}(K), en un conjunto {m_{i,K}(τ)}_{i=1}^{n_{k}^{K}}. Entonces como es usual, sea

{ψ_{j,K}}_{j=1}^{n_{k}^{K}} la base canónica de H_{k}^{K}, tal que, dado i ∈ {1, 2, …, n_{k}^{K}}, se cumple que ψ_{i,K} es el único elemento en H_{k}^{K} tal que:

m_{j,K}(ψ_{i,K}) = δ_{ij}  ∀ j ∈ {1, 2, …, n_{k}^{K}} ,  (18)

donde, en particular, se satisface que:

τ = ∑_{j=1}^{n_{k}^{K}} m_{j,K}(τ) ψ_{j,K} ∀ τ ∈ H_{k}^{K} .

Ahora, para cada K ∈ T_{h} se define la forma bilineal s^{K} : H_{k}^{K} × H_{k}^{K} → ℝ asociada con la matriz identidad en ℝ^{n_{k}^{K}×n_{k}^{K}}, respecto a la base {ψ_{j,K}}_{j=1}^{n_{k}^{K}} de H_{k}^{K}. Esto es:

s^{K}(τ, ζ) := ∑_{i=1}^{n_{k}^{K}} m_{i,K}(τ) m_{i,K}(ζ) ∀ τ, ζ ∈ H_{k}^{K} .

Es gracias a esto que se define S^{K} : ℍ_{h}^{K} × ℍ_{h}^{K} → ℝ como la forma bilineal asociada con los grados de libertad de ℍ_{h}^{K}, tal y como sigue:

S^{K}(τ , ζ) := ∑_{i=1}^{2} s^{K} ((τ_{i1}, τ_{i2})^{t} , (ζ_{i1}, ζ_{i2})^{t}) ,  (19)

para todo τ := (τ_{ij}), ζ := (ζ_{ij}) ∈ ℍ_{h}^{K}.

3.2 Dos elecciones para el proyector Π̂_{k}^{K}

Considere P_{k}^{K} : L^{2}(K) → P_{k}(K) como el proyector ortogonal, para el cual, dado v ∈ L^{2}(K), se caracteriza por:

P_{k}^{K}(v) ∈ P_{k}(K)  y  ∫_{K} P_{k}^{K}(v) q = ∫_{K} v q  ∀ q ∈ P_{k}(K). (20)

Asimismo, sea P_{k}^{K} : L^{2}(K) → P_{k}(K) su correspondiente versión vectorial, tal que para v ∈ L^{2}(K), se cumple que:

P_{k}^{K}(v) ∈ P_{k}(K)  y  ∫_{K} P_{k}^{K}(v)·q = ∫_{K} v·q  ∀ q ∈ P_{k}(K) , (21)

donde, P_{k}^{K}(v) = (P_{k}^{K}(v_{1}), P_{k}^{K}(v_{2}))^{t} para todo v := (v_{1}, v_{2})^{t} ∈ L^{2}(K).

Para efectos de procedimientos posteriores, de acuerdo con [10, Sección 3.2], se tiene que los grados de libertad dados en (15) permiten el cálculo explícito de P_{k}^{K}(τ) para cada τ ∈ H_{k}^{K}. Es decir, es posible determinar el proyector ortogonal para elementos del espacio virtual. En efecto, es suficiente verificar que el lado derecho de (21) se puede calcular en estos casos. Para ello, nótese que gracias a las definiciones de m_{q,e}^{n}(τ) y m_{q,K}^{div}(τ) (ver (15)), es posible determinar el valor de div(τ) ∈ P_{k}(K) empleando la identidad:

∫_{K} div(τ) q = − ∫_{K} τ · ∇q + ∫_{∂K} τ · n q  ∀ q ∈ P_{k}(K) .  (22)

Luego, dado q∈P_{k}(K), se sabe que existen únicos q^{⊥}∈(∇P_{k+1}(K))^{⊥}∩P_{k}(K) y q̃ ∈ P_{k+1}(K), tales que: q = q^{⊥} + ∇q̃. En este sentido, se sigue que:

∫_{K} τ · q = ∫_{K} τ · q^{⊥} + ∫_{K} τ · ∇q̃ = ∫_{K} τ · q^{⊥} − ∫_{K} q̃div(τ) + ∫_{∂K} τ · n q̃,

lo que, en concordancia con (22) y la definición de m_{q,K}^{rot}(τ) (ver (15)), concede el cálculo requerido.

Ahora, de acuerdo con [18, Sección 4.4], es posible considerar al menos dos definiciones para el proyector Π̂_{k}^{K} : ℍ(div; K) → ℍ̂_{k}^{K}, donde ℍ̂_{k}^{K} corresponde a un subespacio polinomial de ℍ(div; K), para cada K ∈ T_{h}. Así, a continuación se definen las dos posibles elecciones para ℍ̂_{k}^{K}, así como la proyección Π̂_{k}^{K}.

3.2.1 Primera elección

Considere primero ℍ̂_{k}^{K} := ℙ_{k}(K), y defina simplemente a Π̂_{k}^{K} como el proyector 𝕃^{2}(K)-ortogonal P_{k}^{K} : 𝕃^{2}(K) → ℙ_{k}(K). En otras palabras, P_{k}^{K} corresponde al operador P_{k}^{K} (ver (21)) actuando en cada fila de un tensor de 𝕃^{2}(K). A su vez, nótese que en el presente caso 2D se cumple que:

dim ℍ̂_{k}^{K} = 4 dim P_{k}(K) = 2(k + 1) (k + 2) .

3.2.2 Segunda elección

En virtud de la identidad σ^{d} = μ ∇u, con div(u) = 0 en Ω (ver (2)), sugiere considerar el proyector definido en [17, Sección 4] para el problema lineal de Stokes. Más precisamente, considere el subespacio de ℙ_{k}(K) dado por:

ℍ̂_{k}^{K} := ℍ̂_{k,∇}^{K} ⊕ ℍ̂_{k,I}^{K} ,  (23)

donde

ℍ̂_{k,∇}^{K} := {∇curl(q) : q ∈ span{x^{α} : 2 ≤ |α| ≤ k + 2} ⊆ P_{k+2}(K)}

y

ℍ̂_{k,I}^{K} := {q I : q ∈ P_{k}(K)} .

Entonces, de [17, Lemma 4.1], es importante recordar que:

dim ℍ̂_{k}^{K} = (k + 1) (k + 4) .

Más aún, observe de (23) que para todo τ ∈ ℍ̂_{k,∇}^{K} se cumple que tr(τ) = 0.

Por otro lado, tal y como se describe en [17], el operador de proyección correspondiente Π̂_{k}^{K} : ℍ(div; K) → ℍ̂_{k}^{K} está definido en términos de la descomposición:

Π̂_{k}^{K}(ζ) := ζ̂_{∇} + q_{ζ} I + c_{ζ} I ∈ ℍ̂_{k}^{K} ,

donde las componentes ζ̂_{∇} ∈ ℍ̂_{k,∇}^{K}, q_{ζ} ∈ P̂_{k}(K) := span{x^{α} : 1 ≤ |α| ≤ k}, y c_{ζ} ∈ ℝ se calculan, secuencialmente, de los siguientes problemas:

• hallar ζ̂_{∇} ∈ ℍ̂_{k,∇}^{K} tal que:

∫_{K} ζ̂_{∇} : τ = ∫_{K} ζ : τ  ∀ τ ∈ ℍ̂_{k,∇}^{K} ,  (24a)

• hallar q_{ζ} ∈ P̂_{k}(K) tal que:

∫_{K} div(q_{ζ} I)·div(q I) = ∫_{K} div(ζ −ζ̂_{∇})·div(q I) ∀ q ∈ P̂_{k}(K) , (24b)

• hallar c_{ζ} ∈ ℝ tal que:

∫_{K} tr(Π̂_{k}^{K}(ζ)) = ∫_{K} tr(ζ) .  (24c)

Respecto a la definición anterior de Π̂_{k}^{K}, observe que la solvencia única de (24) está garantizada en [17, Sección 4]. En particular, nótese que una vez calculado ζ̂_{∇} ∈ ℍ̂_{k,∇}^{K} y luego q_{ζ} ∈ P̂_{k}(K), la identidad (24c) establece que:

c_{ζ} = \frac{1}{2|K|} ∫_{K} \{tr(ζ) − 2q_{ζ}\} .

Ahora, es importante comprobar que los lados derechos de (24) en efecto se pueden calcular para elementos ζ que pertenecen al espacio virtual ℍ_{h}^{K} (ver (9)). Primero que todo, para cada τ := ∇curl(q) ∈ ℍ̂_{k,∇}^{K} se tiene que:

∫_{K} ζ : τ = ∫_{K} ζ : ∇curl(q) = − ∫_{K} curl(q) · div(ζ) + ∫_{∂K} ζn · curl(q) , lo que establece que el lado derecho de (24a) se puede calcular explícitamente para ζ ∈ ℍ_{h}^{K}. Luego, dado que div(ζ) ∈ P_{k}(K) (ver (9)) y ζ̂_{∇} ∈ ℙ_{k}(K), se obtiene que el lado derecho de (24b) también es calculable para cada q ∈ P̂_{k}(K). Finalmente, para el lado derecho de (24c) se sigue que:

∫_{K} tr(ζ) = ∫_{K} ζ : I = ∫_{K} ζ : ∇x = − ∫_{K} x · div(ζ) + ∫_{∂K} ζn · x,

lo cual, de acuerdo con (9), también es calculable. Nótese que el cálculo de Π̂_{k}^{K}(ζ) no requiere de los grados de libertad m_{ρ,K}^{rot}.

4 Ensamblaje del sistema lineal

En esta sección se describe el sistema lineal obtenido a través del problema discreto definido en (10). Para ello, es importante tener en cuenta que la descomposición T_{h} requiere una estructura computacional un poco más sofisticada que las empleadas en los métodos de elementos finitos tradicionales. En efecto, para el caso de los elementos virtuales, y en particular para el problema estudiado en este artículo, lo primero a considerar es que los elementos que conforman T_{h} no necesariamente poseen la misma geometría. Es decir, T_{h} puede estar conformada por elementos de diferentes formas que afectan aspectos importantes, como por ejemplo: la cantidad de aristas que posee, la cantidad de elementos vecinos con los que comparten una aristas, el cálculo de su diámetro, por mencionar algunos.

En particular, cualquier estructura que se implemente para el manejo de la malla, desde un punto de vista de conectividad (topología), debe ser capaz de indicar:

• la cantidad de nodos (puntos),

• la cantidad de aristas en total,

• la cantidad de elementos,

• la cantidad de aristas de frontera,

• las coordenadas de un nodo específico,

• en el caso de un elemento:

– el número de nodos del elemento,

– el índice (identificador) global de un nodo del elemento,

– los índices locales de los nodos en una arista específica, – los índices globales de los nodos en una arista específica,

– el índice global de una arista del elemento,

• en el caso de una arista:

– el índice global de un nodo de la arista,

– el índice global de la arista,

– el índice global del vecino compartido por esta arista,

– los índices globales de los elementos que contienen la arista,

– un identificador que establezca si la arista es de frontera. Asimismo, desde un punto de vista geométrico, para cada K ∈ T_{h}, se debe poder calcular:

• el número de vértices, o equivalentemente, el número de aristas,

• las coordenadas de un vértice específico,

• las coordenadas del baricentro de K,

• el diámetro de K,

• los puntos medios, vectores normales y longitudes de las aristas.

Con ayuda de la información previa, a continuación se detalla la estructura del sistema lineal inducido por (10), así como la forma de ensamblar este. En efecto, empleando (11), (12) y (6) en (10), se obtiene que:

∑_{K∈T_{h}} [\frac{1}{μ} ∫_{K} (Π̂_{k}^{K}(σ_{h}))^{d} : (Π̂_{k}^{K}(τ_{h}))^{d}

+ \frac{1}{α} ∫_{K} div(σ_{h}) · div(τ_{h}) + S^{K} (σ_{h} − Π̂_{k}^{K}(σ_{h}), τ_{h} − Π̂_{k}^{K}(τ_{h})) ]

= ∑_{K∈T_{h}} [−\frac{1}{α} ∫_{K} f · div(τ_{h}) + ∑_{e∈∂K∩Γ} ∫_{e} τ_{h}n · g] ∀ τ_{h} ∈ ℍ_{h} ,

donde Π̂_{k}^{K} : ℍ(div; K) → ℍ̂_{k}^{K} es cualquiera de los dos proyectores definidos en la Sección 3.2. Ahora, con el fin de relajar las restricciones de ℍ_{h} se introduce ξ ∈ ℝ como el multiplicador de Lagrange que permite extraer de ℍ_{h} la condición de integral de traza nula de σ_{h}. Más precisamente, ξ = 0 es una variable

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(2): 215–251, Jul–Dec 2019

auxiliar que permite reescribir el sistema previo de la forma:

∑_{K∈T_{h}} [\frac{1}{μ} ∫_{K}(Π̂_{k}^{K}(σ_{h}))^{d} : (Π̂_{k}^{K}(τ_{h}))^{d} + \frac{1}{α} ∫_{K} div(σ_{h}) · div(τ_{h})

+ S^{K}(σ_{h} − Π̂_{k}^{K}(σ_{h}), τ_{h} − Π̂_{k}^{K}(τ_{h}))

+ ξ ∫_{K} tr(τ_{h})] = ∑_{K∈T_{h}} [−\frac{1}{α} ∫_{K} f · div(τ_{h}) + ∑_{e∈∂K∩Γ} ∫_{e} τ_{h}n · g]

λ ∑_{K∈T_{h}} ∫_{K} tr(σ_{h}) = 0,

para todo (τ_{h}, λ) ∈ ℍ̃_{h} × ℝ, donde ℍ̃_{h} se define como:

ℍ̃_{h} := {τ ∈ ℍ(div; Ω) : τ |_{K} ∈ ℍ_{k}^{K} ∀ K ∈ T_{h}} ,

es decir, ℍ̃_{h} es idéntico a ℍ_{h} (ver (7)), salvo que ya no se impone la condición de integral de traza nula sobre sus elementos.

De acuerdo a lo anterior, el sistema lineal global asociado a (10) posee la estructura:

( \begin{matrix}A + B + H & C \\ C^{t} & 0 \end{matrix} ) ( \begin{matrix}σ & h \\ ξ \end{matrix} ) = ( \begin{matrix}b_{1} + b_{2} \\ 0 \end{matrix} ) ,  (25)

donde:

A:=∑_{K∈T_{h}} \frac{1}{μ}∫_{K}(Π̂_{k}^{K}(σ_{h}))^{d}:(Π̂_{k}^{K}(τ_{h}))^{d},  B:= ∑_{K∈T_{h}}\frac{1}{α}∫_{K}div(σ_{h})·div(τ_{h}),

H:=∑_{K∈T_{h}}S^{K}(σ_{h}−Π̂_{k}^{K}(σ_{h}), τ_{h}−Π̂_{k}^{K}(τ_{h})), C:= ∑_{K∈T_{h}} ξ ∫_{K}tr(τ_{h}),

b_{1} :=∑_{K∈T_{h}} −\frac{1}{α}∫_{K} f · div(τ_{h}),  b_{2} := ∑_{K∈T_{h}} ∑_{e∈∂K∩Γ} ∫_{e}τ_{h}n · g.

Ahora, como es usual en los métodos de elementos finitos, la construcción explícita de la matriz de coeficientes y el vector de la derecha en el sistema (25), se efectúa mediante el ensamblaje de operadores discretos locales provenientes de cada elemento K ∈ T_{h}. Más precisamente, para cada K ∈ T_{h}, considere la versión local del sistema (25) dada por:

( \begin{matrix}A^{K} + B^{K} + H^{K} & C^{K} \\ (C^{K})^{t} & 0 \end{matrix} ) ( \begin{matrix}σ_{h}|_{K} \\ ξ|_{K} \end{matrix} ) = ( \begin{matrix}b_{1}^{K} + b_{2}^{K} \\ 0 \end{matrix} ) ,

donde, recordando de la Sección 3.1 que {ψ_{j,K}}_{j=1}^{n_{k}^{K}} denota la base canónica de

H_{k}^{K}, se define {Ψ_{j,K}}_{j=1}^{2n_{k}^{K}} como la base canónica de ℍ_{k}^{K}, donde

Ψ_{1,K} := ( \begin{matrix}ψ & 1,K \\ & 0 \end{matrix} ) ,  Ψ_{2,K} := ( \begin{matrix}ψ & 2,K \\ & 0 \end{matrix} ) , …,

Ψ_{n_{k}^{K},K} := ( \begin{matrix}ψ & nK_{k},K \\ & 0 \end{matrix} ) ,  Ψ_{n_{k}^{K}+1,K} := ( \begin{matrix}& 0 \\ ψ & 1,K \end{matrix} ) ,

Ψ_{n_{k}^{K}+2,K} := ( \begin{matrix}& 0 \\ ψ & 2,K \end{matrix} ) , …,  Ψ_{2n_{k}^{K},K} := ( \begin{matrix}& 0 \\ ψ & nK_{k},K \end{matrix} ) .

Gracias a esta, se tiene que:

σ_{h}|_{K} = ∑_{j=1}^{2n_{k}^{K}} σ_{h,j}^{K} Ψ_{j,K} ,  (26)

donde {σ_{h,j}^{K}}_{j=1}^{2n_{k}^{K}} son los valores de los grados de libertad de σ_{h}|_{K} ∈ ℍ_{k}^{K} (ver (17)), los cuales a su vez son las incógnitas a determinar. Seguidamente, no es difícil notar que los operadores discretos locales vienen dados por:

• A^{K} := [\frac{1}{μ} ∫_{K} (Π̂_{k}^{K}(Ψ_{i,K}))^{d} : (Π̂_{k}^{K}(Ψ_{j,K}))^{d}] ∈ ℝ^{2n_{k}^{K}×2n_{k}^{K}},

• B^{K} := [\frac{1}{α} ∫_{K} div (Π̂_{k}^{K}(Ψ_{i,K})) · div (Π̂_{k}^{K}(Ψ_{j,K}))] ∈ ℝ^{2n_{k}^{K}×2n_{k}^{K}},

• H^{K} := [S^{K} (Ψ_{i,K} − Π̂_{k}^{K}(Ψ_{i,K}), Ψ_{j,K} − Π̂_{k}^{K}(Ψ_{j,K}))] ∈ ℝ^{2n_{k}^{K}×2n_{k}^{K}},

• C^{K} := [∫_{K} tr(Ψ_{i,K})] ∈ ℝ^{2n_{k}^{K}×1},

• b_{1}^{K} := [−\frac{1}{α} ∫_{K} f · div(Ψ_{i,K})] ∈ ℝ^{2n_{k}^{K}×1},

• b_{2}^{K} := [ ∑_{e∈∂K∩Γ} ∫_{e} Ψ_{i,K}n · g] ∈ ℝ^{2n_{k}^{K}×1}.

La construcción explícita de todos estos operadores discretos se detalla a lo largo de la Sección 5. Por el momento, se considera que, para K ∈ T_{h}, ya estos se han calculado, por lo que a continuación se describe el ensamblaje global del sistema lineal (25). En efecto, sea N el tamaño del sistema (25) el cual viene dado por:

N := \underbrace{2(k + 1)}_{m^{n}}(# de aristas en T_{h})

+ ( \underbrace{(k + 1)(k + 2) − 2}_{m^{div}} + \underbrace{k(k + 1)}_{m^{rot}})(# de elementos en T_{h}) +\underbrace{1}_{ξ} ,

por lo que se tienen N incógnitas asociadas a los grados de libertad. La distribución de estas consiste en primero ordenar los momentos de la componente normal m^{n} para todas las aristas en T_{h}. Estos momentos se ordenan por aristas y no por elementos, con el fin de garantizar continuidad en la componente normal de la solución para que esta, en efecto, pertenezca a ℍ(div; Ω).

Luego, se ordenan los momentos interiores m^{div} y m^{rot} asociados a cada elemento, y finalmente la incógnita correspondiente al multiplicador de Lagrange ξ. Por otro lado, de manera local sobre un K ∈ T_{h} particular, los 2n_{k}^{K} +1 grados de libertad se ordenan tal y como se describe en la Figura 1.

Figura 1: Convención en el orden de los grados de libertad locales para cada K ∈ T_{h}.

De acuerdo con la discusión anterior, suponga que A_{h} ∈ ℝ^{N×N} y b_{h} ∈ ℝ^{N×1} denotan la matriz de coeficientes y el vector de la derecha de (25), respectivamente. Así, estos se pueden ensamblar a través del siguiente algoritmo:

1. Definir A_{h} := 0 y b_{h} := 0

2. Para cada K ∈ T_{h} realizar:

3. Construir los operadores discretos: A^{K}, B^{K}, H^{K}, C^{K}, b_{1}^{K} y b_{2}^{K} 4. Definir m_{K} := \frac{1}{2}(k + 1)(k + 2) − 1 + \frac{1}{2}k(k + 1)

5. Definir n_{K} := (k + 1)d_{K} + m_{K}

6. Definir p ∈ ℝ^{n_{K}×1} y q ∈ ℝ^{n_{K}×1}

7. Para cada e = 1 hasta d_{K} realizar:

8. Definir I_{e} como el índice global en T_{h} de la arista local e

9. Definir w := (I_{e} − 1) · 2(k + 1)

10. Para cada i = 1 hasta k + 1 realizar:

11. Definir s := (e − 1)(k + 1) + i

12. Ingresar p_{s} := w + i

13. Ingresar q_{s} := w + (k + 1) + i

14. Fin del Para de i

15. Fin del Para de e

16. Definir I_{K} como el índice global en T_{h} del elemento K

17. Definir w := (# de aristas en T_{h}) · 2(k + 1) + (I_{K} − 1)(2m_{K})

18. Para cada i = 1 hasta m_{K} realizar:

19. Definir s := (k + 1)d_{K} + i

20. Ingresar p_{s} := w + i

21. Ingresar q_{s} := w + m_{K} + i

22. Fin del Para de i

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 26(2): 215–251, Jul–Dec 2019

23. Definir u := ( \begin{matrix}p \\ q \end{matrix} ) ∈ ℝ^{(2n_{k}^{K})×1}

24. A_{h}(u, u) := A_{h}(u, u) + A^{K}

25. A_{h}(u, N) := A_{h}(u, N) + C^{K}

26. A_{h}(N, u) := A_{h}(N, u) + (C^{K})^{t}

27. b_{h}(u) := b_{h}(u) + b_{1}^{K} + b_{2}^{K}

28. Fin del Para de K

Respecto al procedimiento descrito previamente es importante recordar que la construcción de los operadores discretos de la línea 3 se describe en la Sección 5. Por otro lado, las líneas desde la 6 hasta la 23, construyen un vector u cuyo fin consiste en mapear los grados de libertad locales a su correspondiente ubicación en el sistema global (25). De hecho, p mapea la primera fila del tensor σ_{h}, mientras que q mapea la segunda. Finalmente, el ensamblaje propiamente se realiza en las líneas desde la 24 hasta la 27, donde, en particular, la notación A_{h}(i, j) hace referencia al acceso del bloque de A_{h} obtenido a través de las filas con índices almacenados en i y columnas con índices almacenados en j. El programa MATLAB dispone de este tipo de operaciones de manera sencilla.

5 Construcción de los operadores discretos locales

En esta sección se describe en detalle una forma de construir los operadores discretos locales: A^{K}, B^{K}, H^{K}, C^{K}, b_{1}^{K} y b_{2}^{K}, para cada K ∈ T_{h}. Para ello, se sigue la metodología empleada en [31], en la cual, se construye cada operador a través de la construcción previa de operadores auxiliares intrínsecamente relacionado con el método mixto de elementos virtuales, definido por el subespacio H_{k}^{K} (ver (8)). Para iniciar, es importante aclarar que la base {ψ_{j,K}}_{j=1}^{n_{k}^{K}} de H_{k}^{K} se denomina “virtual” dado que esta no se conoce realmente (por ejemplo, no se sabe si es polinomial), si no solo cómo interactúa esta sobre K. Más precisamente, de la definición de H_{k}^{K} (ver (7)), se intuye que para cada ψ_{j,K} se conocen sus componentes normales, su divergencia y su rotacional. Sin embargo, esto no es del todo preciso. En efecto, lo que realmente se conoce de cada ψ_{j,K} son los valores de sus momentos (15), los cuales están dados en (18). Es claro que gracias a estos es posible determinar a través de cálculos usuales sus componentes normales, su divergencia y su rotacional, explícitamente. Pero para efectos de la siguiente implementación computacional, las identidades (18) son las únicas requeridas.

5.1 Preliminares

De acuerdo a la discusión previa, es necesario conocer funciones polinomiales que permitan realizar cálculos de manera clara. Por ello, sean K ∈ T_{h} y k ≥ 0. Luego, considere ϕ_{1}^{e}, ϕ_{2}^{e}, …, ϕ_{k+1}^{e} la base sobre la arista e ∈ ∂K definida en (13). Es decir,

ϕ_{i}^{e}(x) := (\frac{x − x_{e}}{h_{e}} )^{i−1} , para i = 1, 2, …, k + 1 ,

donde x_{e} es el punto medio de e y h_{e} su longitud. Ahora, sean φ_{1}^{K}, φ_{2}^{K}, …, φ_{M}^{K}, con M := \frac{(k+1)(k+2)}{2} la base (14) dada por:

φ_{i}^{K}(x, y) := (\frac{x − x_{K}}{h_{K}} )^{α} (\frac{y − y_{K}}{h_{K}} )^{β} , para i = 1, 2, …, M ,

donde α + β ∈ {0, 1, …, k}, (x_{K}, y_{K}) el baricentro de K y h_{K} el diámetro de K.

En particular, se asume que el orden corresponde a:

i := \frac{(α + β + 1)(α + β + 2)}{2} − β , con α, β ∈ ℕ y α + β ≤ k ,

esto con el único fin de garantizar la jerarquía de la base, lo que en particular esteblece que φ_{1}^{K} es el polinomio constante. Luego, en virtud de esta base, se define la base de P_{k}(K):

φ_{1}^{K} := ( \begin{matrix}φ_{1}^{K} \\ 0 \end{matrix} ) ,  φ_{2}^{K} := ( \begin{matrix}φ_{2}^{K} \\ 0 \end{matrix} ) , …, φ_{M}^{K} := ( \begin{matrix}φ_{M}^{K} \\ 0 \end{matrix} ) ,
φ_{M+1}^{K} := ( \begin{matrix}0 \\ φ_{1}^{K} \end{matrix} ) , φ_{M+2}^{K} := ( \begin{matrix}0 \\ φ_{2}^{K} \end{matrix} ) , …, φ_{2M}^{K} := ( \begin{matrix}0 \\ φ_{M}^{K} \end{matrix} ) ,

así como la base de ℙ_{k}(K):

Φ_{1}^{K} := ( \begin{matrix}φ_{1}^{K} & 0 \\ 0 & 0 \end{matrix} ) , …, Φ_{M}^{K} := ( \begin{matrix}φ_{M}^{K} & 0 \\ 0 & 0 \end{matrix} ) ,
Φ_{M+1}^{K} := ( \begin{matrix}0 & φ_{1}^{K} \\ 0 & 0 \end{matrix} ) , …, Φ_{2M}^{K} := ( \begin{matrix}0 & φ_{M}^{K} \\ 0 & 0 \end{matrix} ) ,
Φ_{2M+1}^{K} := ( \begin{matrix}0 & 0 \\ φ_{1}^{K} & 0 \end{matrix} ) , …, Φ_{3M}^{K} := ( \begin{matrix}0 & 0 \\ φ_{M}^{K} & 0 \end{matrix} ) ,
Φ_{3M+1}^{K} := ( \begin{matrix}0 & 0 \\ 0 & φ_{1}^{K} \end{matrix} ) , …, Φ_{4M}^{K} := ( \begin{matrix}0 & 0 \\ 0 & φ_{M}^{K} \end{matrix} ) .

Es importante aclarar que las bases {ϕ_{i}^{e}}_{i=1}^{k+1} y {φ_{i}^{K}}_{i=1}^{M} son propias del elemento K, es decir, estas deben construirse para cada elemento K ∈ T_{h}. La razón de esto, es que la geometría de cada K no necesariamente es la misma. Luego, es importante observar que en virtud de la propiedad jerárquica de la base {φ_{i}^{K}}, esta puede ser extendida secuencialmente para grado k +1, lo cual será utilizado a lo largo de esta sección.

Por otro lado, algunas matrices que ayudan en la implementación del esquema (10) de elementos virtuales son las siguientes:

• La matriz de masa en la arista:

M_{mass,1} := [∫_{e} ϕ_{i}^{e}ϕ_{j}^{e}] ∈ ℝ^{(k+1)×(k+1)}.

Con el fin de calcular esta matriz, nótese que:

∫_{e} ϕ_{i}^{e}ϕ_{j}^{e} = h_{e} ∫_{0}^{1} ϕ̂_{i}(x)ϕ̂_{j}(x) dx ,

donde ϕ̂_{i}(x) := (x − \frac{1}{2})^{i−1}, para i = 1, 2, …, k + 1. Así, las entradas de la matriz M_{mass,1} se pueden determinar empleando una regla de cuadratura suficientemente precisa sobre [0, 1]. Además, se cumple que:

M_{mass,1} := h_{e} [∫_{0}^{1} ϕ̂_{i}(x)ϕ̂_{j}(x) dx] ∈ ℝ^{(k+1)×(k+1)},

donde la matriz de la derecha puede ser precomputada independientemente de la arista e. En particular, esto es importante dado que se requiere:

M_{mass,1}^{−1} := \frac{1}{h_{e}} [∫_{0}^{1} ϕ̂_{i}(x)ϕ̂_{j}(x) dx]^{−1} .

• La siguiente matriz:

M_{edge} := [∫_{e} ϕ_{i}^{e}φ_{j}^{K}] ∈ ℝ^{(k+1)× \frac{(k+2)(k+3)}{2}} ,  (27)

donde se utiliza la base {φ_{i}^{K}} hasta grado k + 1. Luego, observe que:

∫_{e} ϕ_{i}^{e}φ_{j}^{K} = h_{e} ∫_{0}^{1} ϕ̂_{i}(x) φ_{j}^{K} ((1 − x) v_{1} + x v_{2}) dx ,

donde v_{1} y v_{2} son los vértices de e, tales que el índice global de v_{1} en T_{h} es menor que el de v_{2}. Esto último garantiza seguir una única orientación sobre la misma arista, independiente del elemento a la que pertenezca.

• La matriz de masa en el elemento:

M_{mass,2}^{(k+1)} := [∫_{K} φ_{i}^{K}φ_{j}^{K}] ∈ ℝ^{\frac{(k+2)(k+3)}{2} × \frac{(k+2)(k+3)}{2}} ,

la cual, por el bien de la implementación, se calcula para grado k + 1. Esto dado que la misma será necesaria, así como su submatriz principal M_{mass,2} ∈ ℝ^{M×M}, cuyo rango de índices es [1, M] × [1, M]. En otras palabras, M_{mass,2}^{(k+1)} es la matriz de masa para la base de P_{k+1}(K), mientras que M_{mass,2} es la matriz de masa para la base de P_{k}(K).

Ahora, para calcular las entradas de M_{mass,2}^{(k+1)} no se utiliza una regla de cuadratura sobre K, sino que se emplea el teorema de la divergencia para convertir la integral de área en la suma de integrales de línea. Más precisamente, empleando el teorema de la divergencia de Gauss, nótese que:

∫_{K}(\frac{x − x_{K}}{h_{K}})^{α}(\frac{y − y_{K}}{h_{K}})^{β} = \frac{h_{K}}{α + β + 2}∫_{K}div(\begin{matrix}(\frac{x−x_{K}}{h_{K}} ) & α+1 & (\frac{y−y_{K}}{h_{K}} ) & β \\ (\frac{x−x_{K}}{h_{K}} ) & α & (\frac{y−y_{K}}{h_{K}} ) & β+1 \end{matrix})

=∑_{e∈∂K} \frac{h_{K}}{α + β + 2} ∫_{e} (\frac{x − x_{K}}{h_{K}})^{α} (\frac{y − y_{K}}{h_{K}})^{β}

[ (\frac{x − x_{K}}{h_{K}} ) n_{1}^{e} + (\frac{y − y_{K}}{h_{K}} ) n_{2}^{e}],  (28)

donde n^{e} := (n_{1}^{e}, n_{2}^{e})^{t} es el vector normal unitario exterior a la arista e. Así, dado que todas las entradas de la matriz M_{mass,2}^{(k+1)} tienen la forma ∫_{K} (\frac{x−x_{K}}{h_{K}} )^{α}(\frac{y−y_{K}}{h_{K}} )^{β}, que es suficiente usar (28) para hallarlas.

• Continuando con la base \{φ_{i}^{K}\}i=1^{\frac{(k+2)(k+3)}{2}} de P_{k+1}(K), se considera la matriz de gradiente:

M_{grad} := [∫_{K} ∇φ_{i+1}^{K} · ∇φ_{j+1}^{K}]  (29)
= [∫_{K} ∂_{x}φ_{i+1}^{K}∂_{x}φ_{j+1}^{K} + ∂_{y}φ_{i+1}^{K}∂_{y}φ_{j+1}^{K}] ∈ ℝ^{M̂×M̂},

donde M̂:= \frac{1}{2}(k + 2)(k + 3) − 1. Nótese que nuevamente se extiende la base {φ_{i}^{K}} hasta grado k + 1, pero eliminando el primer elemento constante. Las entradas de M_{grad} se pueden calcular al utilizar dos veces la fórmula (28).

Seguidamente, para definir las siguientes matrices, es importante observar que {∇φ_{i+1}^{K}}M̂ ; i=1 es una base de ∇P_{k+1}(K). De esta forma, el objetivo ahora es construir una base \{q_{i}^{K}\}_{i=1}^{s_{k}} de G_{k}^{⊥}(K), con s_{k} := \frac{1}{2}k(k + 1). Para ello, se requiere hallar, para cada i ∈ {1, 2, …, s_{k}}, las constantes {α_{j}^{(i)}}_{j=1}^{2M} tales que:

q_{i}^{K} := ∑_{j=1}^{2M} α_{j}^{(i)} φ_{j}^{K}  y  ∫_{K} ∇φ_{i+1}^{K} · q_{j}^{K} = 0 , para i = 1, 2, …, M̂.

Equivalentemente, se requiere resolver los sistemas lineales rectangulares locales:

M_{0}A_{0} = 0,  (30)

donde 0 es la matriz nula de ℝ^{M̂×s_{k}}, y las matrices M_{0} ∈ ℝ^{M̂×(2M)} y A_{0} ∈ ℝ^{(2M)×s_{k}} vienen dadas por:

M_{0} := [∫_{K} ∇φ_{i+1}^{K} · φ_{j}^{K}] ,

y

A_{0} := [α_{j}^{(ℓ)}] = ( \begin{matrix}α_{1}^{(1)} & α_{1}^{(2)} & ⋯ & α_{1}^{(s_{k})} \\ α_{2}^{(1)} & α_{2}^{(2)} & ⋯ & α_{2}^{(s_{k})} \\ ⋮ & ⋮ & ... & ⋮ \\ α_{2M}^{(1)} & α_{2M}^{(2)} & ⋯ & α_{2M}^{(s_{k})} \end{matrix} ) ,

para todo i = 1, 2, …, M̂, j = 1, 2, …, 2M y ℓ = 1, 2, …, s_{k}. Respecto a la construcción de la matriz M_{0} esta es similar a la de M_{grad} (ver (29)), al notar que esta se puede ensamblar a nivel de bloques tal y como se muestra:

M_{0} = [ [∫_{K} ∂_{x}φ_{i+1}^{K}φ_{j}^{K}] [∫_{K} ∂_{y}φ_{i+1}^{K}φ_{j}^{K}] ] ,  (31)

donde ambos bloques poseen tamaños M̂ × M. Por otra parte, para determinar A_{0} se debe resolver el sistema (30), para lo cual, se recomienda utilizar la descomposición QR.

Alternativamente, es importante mencionar que los grados de libertad dados por m_{q,K}^{rot} (ver (15)) y m_{ρ,K}^{rot} (ver (17)) también podrían definirse empleando en lugar de G_{k}^{⊥}(K), la base de cualquier espacio polinomial P̃_{k}(K) que cumpla que:

P_{k}(K) = ∇P_{k+1}(K) ⊕ P̃_{k}(K) .

De esta manera, no se requiere aplicar la descomposición QR, y con ello se pueden calcular estos grados de libertad de una manera más eficiente. Para más detalles, se recomienda consultar [11, Sección 2.1], y particularmente [11, eq. (2.10)].

Las matrices definidas previamente son suficientes para implementar el esquema (10) empleando la primera elección del proyector Π̂_{k}^{K} : ℍ(div; K) → ℍ̂_{k}^{K} (ver Sección 3.2). Es decir, empleando el proyector L^{2}(K)-ortogonal. En el caso de la segunda elección, la correspondiente implementación del proyector definido en (24) requiere además las matrices definidas a continuación:

• Extendiendo la base {φ_{i}^{K}} hasta grado k + 2, y luego descartando los elementos de grado ≤ 2 se obtiene {φ_{i+3}^{K}}_{i=1}^{W}, donde W := \frac{(k+3)(k+4)}{2} − 3 = \frac{1}{2}(k + 1)(k + 6). Así, se define la matriz:

M_{curl} := [∫_{K} ∇curl(φ_{i+3}^{K}) : ∇curl(φ_{j+3}^{K})] ∈ ℝ^{W ×W} ,

donde curl(φ_{i}^{K}) := (∂_{y}φ_{i}^{K}, −∂_{x}φ_{i}^{K})^{t}. Las entradas de M_{curl} se pueden calcular al utilizar tres veces la fórmula (28).

• Ahora, considere:

M_{CB} := [∫_{K} curl(φ_{i+3}^{K}) · φ_{j}^{K}] ∈ ℝ^{W ×2M},

la cual se puede calcular por bloques de la forma:

M_{CB} = [ [∫_{K} ∂_{y}φ_{i+3}^{K}φ_{j}^{K}] [− ∫_{K} ∂_{x}φ_{i+3}^{K}φ_{j}^{K}] ] ,

donde ambos bloques poseen tamaño W × M. Observe que el cálculo de M_{CB} es similar al de M_{0} (ver (31)).

• Sea e ∈ ∂K y con ello:

M_{CE}^{e} := [∫_{e} curl(φ_{i+3}) · ϕ_{j}^{e}] ∈ ℝ^{W ×2(k+1)},

donde {ϕ_{j}^{e}}_{j=1}^{2(k+1)} es la versión vectorial de la base {ϕ_{j}^{e}}_{j=1}^{k+1} definida por:

ϕ_{1}^{e} := ( \begin{matrix}ϕ_{1}^{e} \\ 0 \end{matrix} ) , ϕ_{2}^{e} := ( \begin{matrix}ϕ_{2}^{e} \\ 0 \end{matrix} ) , …, ϕ_{k+1}^{e} := ( \begin{matrix}ϕ_{k+1}^{e} \\ 0 \end{matrix} ) ,
ϕ_{k+2}^{e} := ( \begin{matrix}0 \\ ϕ_{1}^{e} \end{matrix} ) , ϕ_{k+3}^{e} := ( \begin{matrix}0 \\ ϕ_{2}^{e} \end{matrix} ) , …, ϕ_{2(k+1)}^{e} := ( \begin{matrix}0 \\ ϕ_{k+1}^{e} \end{matrix} ) .

Luego, se tiene que:

M_{CE}^{e} = [ [∫_{e} ∂_{y}φ_{i+3}^{K}ϕ_{j}^{e}] [− ∫_{e} ∂_{x}φ_{i+3}^{K}ϕ_{j}^{e}] ] ,

donde ambos bloques poseen tamaño W × (k + 1).

• Seguidamente, considere:

M_{GDGC} := [∫_{K} ∇φ_{i+1}^{K} · div (∇curl(φ_{j+3}^{K}))] ∈ ℝ^{(M−1)×W} ,

la cual se puede calcular empleando dos veces la fórmula (28).

• También son necesarios los vectores fila:

c_{x} := [∫_{K} x̃· φ_{j}^{K}] ∈ ℝ^{1×2M} y c_{x}^{e} := [∫_{e} x̃· ϕ_{j}^{e}] ∈ ℝ^{1×2(k+1)},

donde x̃ := (\frac{x−x_{K}}{h_{K}} , \frac{y−y_{K}}{h_{K}} )^{t}. Ambos vectores se pueden calcular por bloques, tal y como sigue:

c_{x} = [ [∫_{K} (\frac{x − x_{K}}{h_{K}} ) φ_{j}^{K}] [∫_{K} (\frac{y − y_{K}}{h_{K}} ) φ_{j}^{K}] ] ,

y

c_{x}^{e} = [ [∫_{e} (\frac{x − x_{K}}{h_{K}} ) ϕ_{j}^{e}] [∫_{e} (\frac{y − y_{K}}{h_{K}} ) ϕ_{j}^{e}] ] .

• Finalmente, se requiere:

M_{GC} := [∫_{K} Φ_{i}^{K} : ∇curl(φ_{j+3}^{K})] ∈ ℝ^{(4M)×W} ,

cuyas entradas se calculan empleando cuatro veces la fórmula (28).

5.2 El operador de proyección P_{0}

Tal y como se discute en [18] y se evidencia en (12), la forma propuesta de trabajar con los elementos virtuales conlleva la utilización de un operador Π̂_{k}^{K} : ℍ(div; K) → ℍ̂_{k}^{K} que proyecte los elementos de ℍ_{k}^{K} sobre un espacio polinomial adecuado. Esto significa que la implementación necesita de una matriz P_{0} ∈ ℝ^{(4M)×(2n_{k}^{K})} que efectúe dicha proyección. En otras palabras, P_{0} es una matriz de cambio de base, donde los elementos de base virtual tensorial {Ψ_{i,K}}_{i=1}^{2n_{k}^{K}} se aproximan a través de combinaciones lineales de los elementos de la base polinomial tensorial {Φ_{i}^{K}}_{i=1}^{4M} por medio P_{0}. Es importante recordar que n_{k}^{K} fue definido en (16). Más aún, en virtud de la Sección 3.2, a continuación se construye la matriz P_{0} para cada una de las posibilidades descritas.

5.2.1 Primera elección

Recuerde que Π̂_{k}^{K} := P_{k}^{K} : 𝕃^{2}(K) → ℙ_{k}(K) es el proyector 𝕃^{2}(K)-ortogonal. Para este, es necesario construir, secuencialmente, las siguientes matrices:

• P_{grad} := (M_{0})^{t} M_{grad}^{−1} ∈ ℝ^{(2M)×M̂}, con

M := \frac{(k + 1)(k + 2)}{2} y M̂:= \frac{(k + 2)(k + 3)}{2} − 1.

• B_{e} := s^{e} (M_{edge})^{t} M_{mass,1}^{−1} ∈ ℝ^{(M̂+1)×(k+1)}, para cada e ∈ ∂K, donde, si I_{v_{1}} y I_{v_{2}} son los índices globales de los vértices que defienen a e, se tiene que:

s^{e} := \{ \begin{matrix}1 & si I_{v_{1}} < I_{v_{2}}, \\ −1 & en otro caso, \end{matrix}  (32)

representa una orientación para e, las cuales permiten que las integrales involucradas en M_{edge} posean el mismo valor, independientemente del elemento. Más aún, se define la matriz B_{e}^{(r,s)} ∈ ℝ^{(s−r+1)×(k+1)} como la submatriz de B_{e} que contiene todas sus columnas, pero solo las filas desde la r hasta la s, con 1 ≤ r < s ≤ M̂ + 1.

• La matriz de la divergencia:

M_{div} := M_{mass,2}^{−1} B̃ ∈ ℝ^{M×n_{k}^{K}},  (33)

con n_{k}^{K} definido en (16) y B̃ ∈ ℝ^{M×n_{k}^{K}} se define a nivel de bloques como:

B̃ := [ B_{e_{1}}^{(1,M)} B_{e_{2}}^{(1,M)} ⋯ B_{e_{d_{k}}}^{(1,M)}  C̃  0_{M×s_{k}} ] ,

donde 0_{M×s_{k}} es la matriz nula en ℝ^{M×s_{k}}, con s_{k} := \frac{1}{2}k(k + 1), y C̃ ∈ ℝ^{M×(M−1)} definida por:

C̃ := [ \begin{matrix}0 & 1×(M−1) \\ −I_{M−1} \end{matrix} ] ,

con 0_{1×(M−1)} el vector fila de ℝ^{1×(M−1)} y I_{M−1} la matriz identidad de ℝ^{(M−1)×(M−1)}.

• La matriz P_{1} ∈ ℝ^{(2M)×n_{k}^{K}} definida por:

P_{1} := −P_{grad} M̃_{mass} M_{div}

+P_{grad} [B_{e_{1}}^{(2,M̂+1)} B_{e_{2}}^{(2,M̂+1)} ⋯ B_{e_{d_{k}}}^{(2,M̂+1)} 0_{M̂×(M−1+s_{k})} ]

+[0_{(2M)×((k+1)d_{K}+M−1)} (I_{Mass} − P_{grad} M_{0}) A_{0} (A_{0}^{t} I_{Mass} A_{0})^{−1} ],

donde M̃_{mass} ∈ ℝ^{M̂×M} es la submatriz de M_{mass,2}^{(k+1)} que considera las filas de la 2 a la M̂ + 1 y las columnas desde la 1 a la M. Más aún, 0_{M̂×(M−1+s_{k})} y 0_{(2M)×((k+1)d_{K}+M−1)} son las matrices nulas de los es-

pacios ℝ^{M̂×(M−1+s_{k})} y ℝ^{(2M)×((k+1)d_{K}+M−1)}, respectivamente. Por otro lado, se tiene que:

I_{Mass} := I_{2} ⊗ M_{mass,2} ∈ ℝ^{(2M)×(2M)} ,  (34)

donde ⊗ corresponde al producto tensorial de Kronecker.

Finalmente, una vez construidas las matrices previas, se ensambla

P̂_{0} := (I_{2} ⊗ M_{mass,2}^{−1}) P_{1} ∈ ℝ^{(2M)×n_{k}^{K}} ,  (35)

y con ella se tiene que la matriz P_{0} ∈ ℝ^{(4M)×(2n_{k}^{K})} viene dada por:

P_{0} := I_{2} ⊗ P̂_{0} .

5.2.2 Segunda elección

Ahora, se considera la proyección definida en (24), donde similar a la sección previa, se requiere calcular algunas matrices. En este caso, cada matriz está asociada a los problemas (24a), (24b) y (24c), respectivamente.

• La matriz P_{a} := M_{curl}^{−1} P_{2} ∈ ℝ^{W ×(2n_{k}^{K})}, con W := \frac{1}{2}(k + 1)(k + 6), donde P_{2} ∈ ℝ^{W ×(2n_{k}^{K})} está definida por:

P_{2} := −M_{CB} I_{div} + ∑_{e∈∂K} M_{CE}^{e} (I_{2} ⊗ N^{e}) .

Aquí, se tiene que:

I_{div} := I_{2} ⊗ M_{div} ∈ ℝ^{(2M)×(2n_{k}^{K})} ,

donde M_{div} fue definida en (33). Más aún, las matrices N^{e} ∈ ℝ^{(k+1)×n_{k}^{K}} corresponden a:

N^{e} := [ 0_{(k+1)×((I_{e}−1)(k+1))} s^{e} M_{mass,1}^{−1} 0_{(k+1)×ω} ] ,

donde ω := (k + 1)d_{K} − I_{e}(k + 1) + M − 1 + s_{k}, y 1 ≤ I_{e} ≤ d_{K} es el índice local de la artista e.

• La matriz P_{b} := (M̃_{grad})^{−1} (M̃_{0} I_{div} − M_{GDGC} P_{a}) ∈ ℝ^{(M−1)×(2n_{k}^{K})}, donde se considera M̃_{grad} ∈ ℝ^{(M−1)×(M−1)} la submatriz de M_{grad} que posee las filas y columnas con índices en [1, M − 1] × [1, M − 1]. Por otro lado, la matriz M̃_{0} ∈ ℝ^{(M−1)×(2M)} es la submatriz de M_{0} que contiene todas sus columnas, pero solo las filas desde la 1 hasta la M − 1.

• La matriz P_{c} ∈ ℝ^{1×(2n_{k}^{K})}, es el vector fila definido por:

P_{c} := \frac{1}{|K|} (−\frac{h_{K}}{2} c_{x} I_{div} − c_{0} P_{b} + \frac{h_{K}}{2} ∑_{e∈∂K} c_{x}^{e} (I_{2} ⊗ N^{e})) ,

donde |K| es el área de K, la cual se encuentra almacenada en la entrada 1-1 de M_{mass,2}. Además, c_{0} ∈ ℝ^{1×(M−1)} es la primera fila de M_{mass,2}, excepto por su primera entrada. Es decir, contiene las columnas desde la 2 hasta la M.

Con ayuda de lo anterior, se concluye la matriz P_{0} ∈ ℝ^{(4M)×(2n_{k}^{K})} de la forma:

P_{0} := (I_{4} ⊗ M_{mass,2}^{−1}) P_{3} ,  (36)

donde P_{3} ∈ ℝ^{(4M)×(2n_{k}^{K})} viene dada por:

P_{3} := M_{GC} P_{a} + [( \begin{matrix}1 \\ 0 \\ 0 \\ 1 \end{matrix} ) ⊗ M̃_{mass}] P_{b} + [( \begin{matrix}1 \\ 0 \\ 0 \\ 1 \end{matrix} ) ⊗ c_{1}] P_{c} ,

donde M̃_{mass} ∈ ℝ^{M̂×M} es la submatriz de M_{mass,2}^{(k+1)} que considera las filas de la 2 a la M̂+ 1 y las columnas desde la 1 a la M, mientras que el vector c_{1} ∈ ℝ^{M×1} es la primera columna de M_{mass,2}.

5.3 El operador A^{K}

El operador A^{K} ∈ ℝ^{(2n_{k}^{K})×(2n_{k}^{K})} se define como la matriz:

A^{K} := [\frac{1}{μ} ∫_{K} (Π̂_{k}^{K}(Ψ_{i,K}))^{d} : (Π̂_{k}^{K}(Ψ_{j,K}))^{d}] = (P_{0})^{t} (M_{μ} ⊗M_{mass,2}) P_{0} ,

donde

M_{μ} := \frac{1}{2μ} ( \begin{matrix}1 & 0 & 0 & −1 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ −1 & 0 & 0 & 1 \end{matrix} ) ,

es la matriz asociada al operador desviador.

5.4 El operador B^{K}

El operador B^{K} ∈ ℝ^{(2n_{k}^{K})×(2n_{k}^{K})} viene dado por:

B^{K} := [\frac{1}{α} ∫_{K} div (Π̂_{k}^{K}(Ψ_{i,K})) · div(Π̂_{k}^{K} (Ψ_{j,K}))]

= \frac{1}{α} I_{2} ⊗ ((M_{div})^{t} M_{mass,2} M_{div}) ,

donde M_{div} ∈ ℝ^{M×n_{k}^{K}} fue definida en (33).

5.5 El operador H^{K}

A diferencia de B^{K}, la construcción del operador H^{K} ∈ ℝ^{(2n_{k}^{K})×(2n_{k}^{K})} definido por:

H^{K} := [S^{K} (Ψ_{i,K} − Π̂_{k}^{K}(Ψ_{i,K}), Ψ_{j,K} − Π̂_{k}^{K}(Ψ_{j,K}))] ,

varía dependiendo del proyector que se emplee (ver Sección 3.2). Por ende, a continuación se discute por separado el ensamblaje de la matriz H^{K} en cada caso.

5.5.1 Primera elección

Para iniciar, en virtud de (19), se cumple que:

H^{K} = I_{2} ⊗ (H_{0}^{t} H_{0}) ,

donde H_{0} ∈ ℝ^{n_{k}^{K}×n_{k}^{K}} se descompone de la forma: H_{0} = H_{1} − H_{2} P̂_{0}, donde P̂_{0} ∈ ℝ^{(2M)×n_{k}^{K}} está definida en (35). Ahora, H_{1} ∈ ℝ^{n_{k}^{K}×n_{k}^{K}} es la matriz diagonal:

H_{1} := ( \begin{matrix}s^{e_{1}} I_{k+1} \\ & s^{e_{2}} I_{k+1} \\ & & . & . & . \\ & & & & & s^{e_{d_{K}}} I_{k+1} \\ 0 & M−1+s_{k} \end{matrix} ) ,  (37)

donde la constante s^{e} se define en (32), mientras que s_{k} := \frac{1}{2}k(k + 1). Por otro lado, considere, para cada e ∈ ∂K, las matrices E^{e} ∈ ℝ^{(k+1)×(2M)} dadas por:

E^{e} := n^{e} ⊗ M̃_{edge} ,

donde n^{e} ∈ ℝ^{1×2} es el vector (fila) normal unitario exterior a la arista e y M̃_{edge} ∈ ℝ^{(k+1)×M} es la submatriz de M_{edge} que contiene todas sus filas, pero las columnas desde la 1 hasta la M.

De esta forma, se concluye que H_{2} ∈ ℝ^{n_{k}^{K}×(2M)} viene dada por:

H_{2} := ( E^{e_{1}} ; E^{e_{2}} ; ⋮ ; E^{e_{d_{K}}} ; 0_{(M−1+s_{k})×(2M)} ) .  (38)

5.5.2 Segunda elección

De manera similar, para el proyector (24), se escribe:

H^{K} = (H_{CG})^{t} H_{CG} ,

donde H_{CG} := (I_{2} ⊗ H_{3}) − (I_{2} ⊗ H_{4}) P_{0}, con P_{0} ∈ ℝ^{(4M)×(2n_{k}^{K})} definida en (36). Luego, la matriz H_{3} ∈ ℝ^{n_{k}^{K}×n_{k}^{K}} es similar a H_{1} (ver (37)) salvo en su último bloque. Más precisamente, se tiene que:

H_{3} := ( \begin{matrix}s^{e_{1}} I_{k+1} \\ & s^{e_{2}} I_{k+1} \\ & & . & . & . \\ & & & & & s^{e_{d_{K}}} I_{k+1} \\ & & & & & & I & M−1+s_{k} \end{matrix} ) .

Asimismo, la matriz H_{4} ∈ ℝ^{n_{k}^{K}×(2M)} sigue la construcción de H_{2} (ver (38)) tal y como se muestra:

H_{4} := ( E^{e_{1}} E^{e_{2}} ⋮ E^{e_{d_{K}}} M̃_{0} (A_{0})^{t} I_{Mass} ) ,

donde M̃_{0} ∈ ℝ^{(M−1)×(2M)} es la submatriz de M_{0} que contiene todas sus columnas, pero las filas desde la 1 hasta la M − 1. La matriz I_{Mass} ∈ ℝ^{(2M)×(2M)} fue definida en (34).

5.6 El operador C^{K}

El operador C^{K} ∈ ℝ^{(2n_{k}^{K})×1} está definido por:

C^{K} := [∫_{K} tr(Ψ_{i,K})] = (P_{0})^{t} [( \begin{matrix}1 \\ 0 \\ 0 \\ 1 \end{matrix} ) ⊗ c_{1}] ,

donde el vector c_{1} ∈ ℝ^{M×1} es la primera columna de M_{mass,2}.

5.7 El operador b_{1}^{K}

Considere f := (f_{1}, f_{2})^{t} la función fuente, y con ella se calculan:

P_{k}(f_{1}) = ∑_{i=1}^{M} a_{i}^{f_{1}}φ_{i}^{K}  y  P_{k}(f_{2}) = ∑_{i=1}^{M} a_{i}^{f_{2}}φ_{i}^{K} ,  (39)

donde P_{k}^{K} : L^{2}(K) → P_{k}(K) es el proyector L^{2}(K)-ortogonal (ver (20)). Para ello, se emplea la propiedad (20) con la base {φ_{i}^{K}}_{i=1}^{M}, donde es requerida una cuadratura adecuada para dominios poligonales. Algunas de estas se pueden encontrar, por ejemplo, en [32], [30] y [21]. De esta forma, si a^{f_{1}} ∈ ℝ^{M×1} y a^{f_{2}} ∈ ℝ^{M×1} son los vectores coordenadas de f_{1} y f_{2} (que satisfacen (39)), respectivamente, se cumple que el operador b_{1}^{K} ∈ ℝ^{(2n_{k}^{K})×1} viene dado por:

b_{1}^{K} := [−\frac{1}{α} ∫_{K} f · div(Ψ_{i,K})] = −\frac{1}{α} [ \begin{matrix}(M_{div})^{t} a^{f_{1}} \\ (M_{div})^{t} a^{f_{2}} \end{matrix} ] ,

donde M_{div} ∈ ℝ^{M×n_{k}^{K}} está definida en (33).

5.8 El operador b_{2}^{K}

Considere g := (g_{1}, g_{2})^{t} el dato de Dirichlet y defina, para cada e ∈ ∂K, los vectores:

b_{g_{ℓ}}^{e} := \{ \begin{matrix}[∫_{e} g_{ℓ}ϕ_{i}^{e}] = [h_{e} ∫_{0}^{1} g_{ℓ} ((1 − x) v_{1} + x v_{2}) ϕ̂_{i}(x) dx] & si e ∈ Γ , \\ 0_{(k+1)×1} & si e ∉ Γ , \end{matrix}

para ℓ = 1, 2, los cuales son de tamaño (k+1)×1, donde v_{1} y v_{2} son los vértices de e, tales que el índice global de v_{1} en T_{h} es menor que el de v_{2}. Nótese que la construcción de estos vectores es similar a la de las matrices M_{edge} (ver (27)). Ahora, el vector b_{2}^{K} ∈ ℝ^{(2n_{k}^{K})×1} viene dado por:

b_{2}^{K} := [ ∑_{e∈∂K∩Γ} ∫_{e} Ψ_{i,K}n · g] = [ b_{g_{1}}^{e_{1}} ; b_{g_{1}}^{e_{2}} ; ⋮ ; b_{g_{1}}^{e_{d_{K}}} ; 0_{(M−1+s_{k})×1} ; b_{g_{2}}^{e_{1}} ; b_{g_{2}}^{e_{2}} ; ⋮ ; b_{g_{2}}^{e_{d_{K}}} ; 0_{(M−1+s_{k})×1} ] .

5.9 Aproximaciones calculables de σ, u y p

En virtud de los operadores descritos en las secciones previas, es posible seguir la Sección 4 para ensamblar el sistema lineal global, (25) que permite determinar los grados de libertad de σ_{h}. Es importante recordar que estos siguen un ordenamiento donde, primero se encuentran los momentos de la componente normal m^{n} para todas las aristas en T_{h}, luego se ordenan los momentos interiores m^{div} y m^{rot} asociados a cada elemento, respectivamente, y finalmente la incógnita correspondiente al multiplicador de Lagrange ξ. Así, una vez resulto el sistema lineal (25) se pueden recuperar, para cada K ∈ T_{h}, los grados de libertad locales s^{K} := (σ_{h,1}^{K}, σ_{h,2}^{K}, …, σ_{h,n_{k}^{K}}^{K})^{t} (ver (26)), ordenados tal y como se indica en la Figura 1.

De acuerdo a lo anterior, para cada K ∈ T_{h}, es posible considerar la aproximación local de σ ∈ ℍ dada por:

σ̂_{h}^{K} := Π̂_{k}^{K}(σ_{h}|_{K}) ,

la cual es conocida completamente en todo K. Más precisamente, se tiene que:

σ̂_{h}^{K} = ∑_{i=1}^{4M} α_{i}^{K}Φ_{i}^{K} ,  (40)

donde α^{K} := (α_{1}^{K}, α_{2}^{K}, …, α_{4M}^{K})^{t} se obtiene de realizar el producto:

α^{K} = P_{0} s^{K} ,

donde P_{0} ∈ ℝ^{(4M)×(2n_{k}^{K})} está definida en la Sección 5.2.

Finalmente, empleando (40), (4) y (3), se pueden determinar, para cada K ∈ T_{h}, aproximaciones calculables para u y p, dadas por:

u_{h}^{K} := \frac{1}{α}\{P_{k}^{K}(f) + div(σ_{h}|_{K})\}  y  p_{h}^{K} := −\frac{1}{2} tr(σ̂_{h}^{K}) ,

tal y como se introduce en [18, Sección 5.2], donde P_{k}^{K} : L^{2}(K) → P_{k}(K) está definido en (21).

6 Conclusiones y direcciones futuras

En este artículo se ha descrito en detalle una implementación computacional para el método mixto de elementos virtuales presentado en [18] para el problema lineal de Brinkman en dos dimensiones con condiciones de Dirichlet no homogéneas. El trabajo realizado en [18] conlleva un gran detalle en el análisis teórico del método, pero no incluye mayores especificaciones sobre cómo se realizó su programación en MATLAB, a pesar de mostrar experimentos numéricos. Así, el trabajo desarrollado en este escrito complementa [18] en un área bastante importante de detallar, dado que al no contar con un elemento de referencia que agilice las integraciones numéricas, además, del hecho de no conocer explícitamente la base canónica virtual, que la implementación de estos métodos puede considerarse de gran dificultad a primera vista. Sin embargo, la misma consiste en construir operadores discretos adecuados y realizar un ensamblaje global como es usual en los métodos de elementos finitos.

Por otro lado, es importante notar que la cantidad de matrices locales necesarias para el ensamblaje global (ver Sección 5), no necesariamente implica que la implementación propuesta aquí sea poco eficiente. Esto dado que los tamaños de estas matrices, las cuales dependen del grado polinomial k ≥ 0 empleado, son considerablemente pequeños. Más aún, muchas de estas matrices no se construyen si se emplean un bajo orden de aproximación. En este sentido, es importante resaltar que implementaciones de alto orden para los métodos de elementos finitos no son frecuentes de realizar, dado que las mismas se consideran complicadas, sin embargo, aquí se mostró que no es realmente un inconveniente significativo.

Respecto a la validación del código, esta se presenta en la Sección 6 de [18], donde se demuestra que los resultados coinciden con el análisis desarrollado a lo largo de [18]. Además, en [20] se sigue el mismo paradigma de programación para el problema de elasticidad lineal, mientras que en [19] se extiende a un problema no lineal. En este sentido, en [26] se emplearon todos los operadores discretos presentados en la Sección 5 para una versión no lineal del problema de Brinkman estudiado en este trabajo.

Para finalizar, en relación a los trabajos futuros, en particular es de nuestro interés detallar la implementación del esquema virtual introducido en [25] para las ecuaciones de Navier-Stokes, donde se emplean, simultáneamente, subespacios virtuales H^{1}-conformes, junto con subespacios virtuales ℍ(div)-conformes. Es decir, además de describir el método no lineal, se introducen dos subespacios virtuales que permiten implementar, por el momento, la mayoría de problemas de interés estudiados en la literatura sobre métodos de elementos virtuales. Por otro lado, también es de interés extender lo realizado aquí para problemas en tres dimensiones.

Agradecimientos

Este trabajo fue financiado por la Universidad Nacional, a través del proyecto 0106-16.

Referencias

[1] V. Anaya, G. N. Gatica, D. Mora, R. Ruiz-Baier„ An augmented velocityvorticity-pressure formulation for the Brinkman equations, Internat. J. Numer. Methods Fluids 79 (2015), no. 3, 109–137.

[2] V. Anaya, D. Mora, C. Reales, R. Ruiz-Baier, Stabilized mixed approximation of axisymmetric Brinkman flow, ESAIM Math. Model. Numer. Anal. 49 (2015), no. 3, 855–874.

[3] V. Anaya, D. Mora, R. Oyarzúa, R. Ruiz-Baier, A priori and a posteriori error analysis of a mixed scheme for the Brinkman problem, Numer. Math. 133 (2016), no. 4, 781–817.

[4] P.F. Antonietti, L. Beirão da Veiga, D. Mora, M. Verani, A stream virtual element formulation of the Stokes problem on polygonal meshes, SIAM J. Numer. Anal. 52 (2014), no. 1, 386–404.

[5] E. Artioli, L. Beirão da Veiga, C. Lovadina, E. Sacco, E.Arbitrary order 2D virtual elements for polygonal meshes: part I, elastic problem, Comput. Mech. 60 (2017), no. 3, 355–377.

[6] L. Beirão da Veiga, F. Brezzi, A. Cangiani, L. D. Marini, G. Manzini, A. Russo, Basic principles of virtual elements methods, Math. Models Methods Appl. Sci. 23 (2013), no. 1, 199–214.

[7] L. Beirão da Veiga, F. Brezzi, L.D. Marini, Virtual elements for linear elasticity problems, SIAM J. Numer. Anal. 51 (2013), no. 2, 794–812.

[8] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, The hitchhiker’s guide to the virtual element method, Math. Models Methods Appl. Sci. 24 (2014), no. 8, 1541–1573.

[9] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, H(div) and H(curl)-conforming virtual element methods, Numer. Math. 133 (2016), no. 2, 303–332.

[10] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM Math. Model. Numer. Anal. 50 (2016), no. 3, 727–747.

[11] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo,Serendipity face and edge VEM spaces, arXiv preprint, arXiv: 1606.01048v1, 2016.

[12] L. Beirão da Veiga, C. Lovadina, G. Vacca, Virtual elements for the Navier-Stokes problem on polygonal meshes, SIAM J. Numer. Anal. 56 (2018), no. 3, 1210–1242.

[13] F. Brezzi, A. Cangiani, G. Manzini, A. Russo, Mimetic finite differences and virtual element methods for diffusion problems on polygonal meshes, I.M.A.T.I.-C.N.R, Technical Report 22PV12/0/0, (2012), 1–27.

[14] F. Brezzi, R.S. Falk, L.D. Marini, Basic principles of mixed virtual element methods, ESAIM Math. Model. Numer. Anal. 48 (2014), no. 4, 1227–1240.

[15] F. Brezzi, L.D. Marini, Virtual element methods for plate bending problems, Comput. Methods Appl. Mech. Engrg. 253 (2013), 455–462.

[16] E. Cáceres, Mixed virtual element methods. Applications in fluid mechanics, tesis, Ingeniero Civil Matemático, Universidad de Concepción, 2015.

[17] E. Cáceres, G.N. Gatica, A mixed virtual element method for the pseudostress-velocity formulation of the Stokes problem, IMA J. Numer. Anal. 37 (2017), no. 1, 296–331.

[18] E. Cáceres, G.N. Gatica, F.A. Sequeira, A mixed virtual element method for the Brinkman problem, Math. Models Methods Appl. Sci. 27 (2017), no. 4, 707–743.

[19] E. Cáceres, G.N. Gatica, F.A. Sequeira, A mixed virtual element method for quasi-Newtonian Stokes flows, SIAM J. Numer. Anal. 56 (2018), no. 1, 317–343.

[20] E. Cáceres, G.N. Gatica, F.A. Sequeira, A mixed virtual element method for a pseudostress-based formulation of linear elasticity, Applied Numerical Mathematics 135 (2019), 423–442.

[21] E.B. Chin, J.B. Lasserre, N. Sukumar, Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra, Comput. Mech. 56 (2015), no. 6, 967–981.

[22] A.L. Gain, C. Talischi, G.H. Paulino, On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes, Comput. Methods Appl. Mech. Engrg. 282 (2014), 132–160.

[23] G.N. Gatica, L.F. Gatica, A. Márquez, Analysis of a pseudostress-based mixed finite element method for the Brinkman model of porous media flow, Numer. Math. 126 (2014), no. 4, 635–677.

[24] G.N. Gatica, L.F. Gatica, F.A. Sequeira, Analysis of an augmented pseudostress-based mixed formulation for a nonlinear Brinkman model of porous media flow, Comput. Methods Appl. Mech. Engrg. 289 (2015), 104–130.

[25] G.N. Gatica, M. Munar, F.A. Sequeira, A mixed virtual element method for the Navier-Stokes equations, Math. Models Methods Appl. Sci. 28 (2018), no. 14, 2719–2762.

[26] G.N. Gatica, M. Munar, F.A. Sequeira, A mixed virtual element method for a nonlinear Brinkman model of porous media flow, Calcolo 55 (2018), no. 2, article 21.

[27] J. Guzmán, M. Neilan, A family of nonconforming elements for the Brinkman problem, IMA J. Numer. Anal. 32 (2012), no. 4, 1484–1508.

[28] M. Juntunen, R. Stenberg, Analysis of finite element methods for the Brinkman problem, Calcolo 47 (2010), no. 3, 129–147.

[29] J. Könnö, R. Stenberg, H(div)-conforming finite elements for the Brinkman problem, Math. Models Methods Appl. Sci. 21 (2011), no.11, 2227–2248.

[30] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons, Comput. Mech. 47 (2011), no. 5, 535–554.

[31] F.A. Sequeira, P.E. Castillo, Implementación del método LDG para mallas no estructuradas en 3D, Revista de Matemática: Teoría y Aplicaciones 19 (2012), no. 2, 141–156.

[32] A. Sommariva, M. Vianello, Product Gauss cubature over polygons based on Green’s integration formula, BIT Numerical Mathematics 47 (2007), no. 2, 441–453.

[33] P. Vassilevski, U. Villa, A mixed formulation for the Brinkman problem, SIAM J. Numer. Anal 52 (2014), no. 1, 258–281.