Reconstructed 1D basis functions for `P2r/P2r`

In this work, the reconstructed basis is not introduced analytically first. Instead, we define a linear reconstruction operator from vertex degrees of freedom to full quadratic degrees of freedom. The reconstructed basis functions are then obtained from the columns of this operator.

Reconstruction procedure

function midpoints_periodic(x::AbstractVector, u::AbstractVector; chi=0.0)
    N = length(x)
    ux, uxx = recover_derivatives_periodic(x, u)
    um = zeros(eltype(u), N)
    for e in 1:N
        ep = modp(e, N) + 1
        h = (e < N) ? (x[e+1] - x[e]) : (1.0 - x[e] + x[1])
        UL = u[e]  + (h/2)*ux[e]  + (h^2/8)*uxx[e]
        UR = u[ep] - (h/2)*ux[ep] + (h^2/8)*uxx[ep]
        um[e] = (0.5 + chi)*UL + (0.5 - chi)*UR
    end
    return um
end

function build_reconstruction_matrix_periodic(x::AbstractVector; chi=0.0)
    N = length(x)
    nv = N
    ndof_full = 2N
    vid(j) = modp(j, N) + 1
    midid(e) = N + e

    rows = Int[]
    cols = Int[]
    vals = Float64[]

    for j in 0:N-1
        push!(rows, vid(j))
        push!(cols, j+1)
        push!(vals, 1.0)
    end

    for k in 1:nv
        ek = zeros(Float64, nv)
        ek[k] = 1.0
        um = midpoints_periodic(x, ek; chi=chi)
        for e in 1:N
            wk = um[e]
            if abs(wk) > 0.0
                push!(rows, midid(e))
                push!(cols, k)
                push!(vals, wk)
            end
        end
    end

    return sparse(rows, cols, vals, ndof_full, nv)
end

1. Periodic 1D setting

Consider a periodic mesh with vertices

x_1, \dots, x_N,

and cells

K_e = [x_e, x_{e+1}],

with periodic indexing. On the last cell, we have

K_N = [x_N, x_1 + 1]

if the periodic domain is identified with $[0,1]$ .

In one space dimension, a continuous quadratic finite element space has:

one degree of freedom at each vertex,
one degree of freedom at each cell midpoint.

Hence, the full quadratic representation has $2N$ degrees of freedom.

2. Reconstruction of midpoint values

Let

U = (U_1, \dots, U_N)^{\top} \in \mathbb{R}^N

denote the vector of vertex values.

From these values, we first recover approximations of the first and second derivatives at the vertices,

u_x(x_i), \qquad u_{xx}(x_i),

through

ux, uxx = recover_derivatives_periodic(x, u)

On each cell $K_e = [x_e, x_{e+1}]$ of size

h_e = x_{e+1} - x_e.

Then, we form two quadratic Taylor approximations of the midpoint value.

The left-based approximation is

U_L^{(e)} = U_e + \frac{h_e}{2} u_x(x_e) + \frac{h_e^2}{8} u_{xx}(x_e),

while the right-based approximation is

U_R^{(e)} = U_{e+1} - \frac{h_e}{2} u_x(x_{e+1}) + \frac{h_e^2}{8} u_{xx}(x_{e+1}).

The reconstructed midpoint value is then defined by

U_{m,e} = \left(\frac{1}{2}+\chi\right) U_L^{(e)} + \left(\frac{1}{2}-\chi\right) U_R^{(e)}.

For $\chi = 0$ , this gives the symmetric average of the two midpoint predictions. For $\chi \neq 0$ , the reconstruction is biased.

3. Reconstruction operator

The previous construction is linear with respect to the vertex vector $U$ . Therefore, it defines a linear operator

R \in \mathbb{R}^{2N \times N}

such that

u = R U,

where $u \in \mathbb{R}^{2N}$ is the vector of full quadratic degrees of freedom,

u = \begin{pmatrix} U_1 \\ \vdots \\ U_N \\ U_{m,1} \\ \vdots \\ U_{m,N} \end{pmatrix}.

The first $N$ rows of $R$ simply copy the vertex values, while the last $N$ rows provide the reconstructed midpoint values. Hence, $R$ can be written as

R = \begin{pmatrix} I \\ W \end{pmatrix},

where:

$I \in \mathbb{R}^{N \times N}$ is the identity matrix,
$W \in \mathbb{R}^{N \times N}$ maps vertex values to reconstructed midpoint values.

Equivalently,

u = \begin{pmatrix} U \\ W U \end{pmatrix}.

4. Construction of the matrix $R$

Since the reconstruction is linear, the matrix $R$ is obtained by applying the reconstruction to the canonical basis vectors of $\mathbb{R}^N$ .

Let

e^{(k)} = (0,\dots,0,1,0,\dots,0)^{\top}

be the $k$ -th canonical basis vector. Then the $k$ -th column of $R$ is given by the reconstructed full quadratic vector associated with $e^{(k)}$ .

In the code, this is done by

ek = zeros(Float64, nv)
ek[k] = 1.0
um = midpoints_periodic(x, ek; chi=chi)

The first block of the column is $e^{(k)}$ itself, since the vertex values are copied. The second block is the vector of reconstructed midpoint values obtained from $e^{(k)}$ .

Thus, the $k$ -th column of $R$ is

R e^{(k)} = \begin{pmatrix} e^{(k)} \\ W e^{(k)} \end{pmatrix}.

This is exactly what the routine build_reconstruction_matrix_periodic assembles.

5. Reconstructed basis functions

Let $\{\varphi_j\}_{j=1}^{2N}$ denote the standard global continuous $P^2$ basis, ordered as follows:

$\varphi_1, \dots, \varphi_N$ are the vertex basis functions,
$\varphi_{N+1}, \dots, \varphi_{2N}$ are the midpoint basis functions.

Any function in the full continuous $P^2$ space can be written as

u_h(x) = \sum_{j=1}^{2N} u_j \varphi_j(x),

where $u \in \mathbb{R}^{2N}$ is the vector of full quadratic degrees of freedom.

In the reconstructed formulation, the dofs are the vertex values $U \in \mathbb{R}^N$ , and the full quadratic coefficients satisfy

u = R U.

Therefore,

u_h(x) = \sum_{j=1}^{2N} (R U)_j \varphi_j(x).

Reordering the sum with respect to the vertex unknowns $U_k$ gives

u_h(x) = \sum_{k=1}^N U_k \psi_k(x),

where the reconstructed basis functions are defined by

\psi_k(x) = \sum_{j=1}^{2N} R_{j,k} \, \varphi_j(x).

This is the definition of the reconstructed P2r basis.

Using the block decomposition $R = \binom{I}{W}$ , we immediately get

\psi_k = \varphi_k + \sum_{e=1}^N W_{e,k} \, \varphi_{N+e}.

Hence, each reconstructed basis function consists of:

its vertex basis contribution $\varphi_k$ ,
plus a linear combination of midpoint basis functions induced by the reconstruction.

6. Reconstructed space

The reconstructed space is the subspace of the full continuous $P^2$ space for which midpoint degrees of freedom are constrained by the reconstruction from vertex values. It can be written as

V_h^r = \left\{ v_h \in V_h^{P^2} \;:\; \text{the midpoint dofs of } v_h \text{ are reconstructed from its vertex dofs} \right\}.

A natural basis of this space is given by the reconstructed basis functions $\{\psi_k\}_{k=1}^N$ .

7. Matrix form of reconstructed operators

Let $A \in \mathbb{R}^{2N \times 2N}$ be any matrix assembled in the full quadratic $P^2$ space, for instance a mass, stiffness, or advection matrix.

Since the reconstructed solution satisfies $u = R U$ , the corresponding matrix in the reconstructed vertex space is

A^r = R^T A R.

Indeed, for two reconstructed vectors $u = R U$ and $v = R V$ , we have

v^T A u = (R V)^T A (R U) = V^T (R^T A R) U.

Thus, the reconstructed operator is obtained by restricting the full operator to the reconstructed subspace. For a mixed operator, we insert $R$ on the side corresponding to the reconstructed space. If the trial space is reconstructed and the test space is unchanged, then

A^r = A R.

If the test space is reconstructed and the trial space is unchanged, then

A^r = R^T A.

If both trial and test spaces are reconstructed, then

A^r = R^T A R.

from a full P2/P2 matrix $A_{P2/P2}$ , we have $A_{P2r/P2r} = R^T A_{P2/P2} R,$
for a mixed operator, we insert $R$ on the trial and/or test side as appropriate, e.g. $A_{P2r/P1} = A_{P2/P1} R.$

Summary

Starting from vertex values $U \in \mathbb{R}^N$ :

recover derivative approximations at the vertices,
use left and right quadratic Taylor expansions to predict midpoint values,
combine the two predictions with parameter $\chi$ ,
obtain a linear reconstruction operator $R : \mathbb{R}^N \to \mathbb{R}^{2N}, \qquad u = R U,$
define the reconstructed basis functions by $\psi_k = \sum_{j=1}^{2N} R_{j,k} \, \varphi_j,$
build reconstructed matrices through $A^r = R^T A R.$

In the code, the reconstructed basis functions are never introduced explicitly. They are the basis functions induced by the reconstruction operator.