Problem initializing a VectorFunctionSpace in Firedrake

[afgalindo]

I am using Firedrake to define a 2-dimensional vector function on a 1D mesh. I want to initialize it with CG1 elements using (x, x) as initial data. Here’s my minimal working example:

from firedrake import *
import numpy as np

Lx = 1.0
Nx = 50  # example number of intervals
mesh = IntervalMesh(Nx, Lx)
x = SpatialCoordinate(mesh)
W = VectorFunctionSpace(mesh, "CG", 1, dim=2)
U = Function(W, name="Solution")

# Attempted initialization
U.interpolate(as_vector((x[0], x[0])))

# Extract components
Er_func = U.sub(0)
F_func = U.sub(1)

# Evaluate at points
x_vals = np.linspace(0, Lx, 100)
        
numerical_Er_vals = np.array([Er_func.at(xi) for xi in x_vals])
numerical_F_vals = np.array([F_func.at(xi) for xi in x_vals])

        # Plot Er(x)
        plt.figure(figsize=(8, 5))
        plt.plot(x_vals, numerical_Er_vals, label="CG $E_r$", linewidth=2)
        plt.xlabel("x")
        plt.ylabel(r"$E_r(x)$")
        plt.title(f"$E_r$ at t = {self.Tmax:.4f}")
        plt.legend()
        plt.grid(True)
        plt.tight_layout()
        plt.savefig(os.path.join(plot_dir, f"Er_vs_exact_t_{self.Tmax:.4f}.png"))
        plt.close()

        # Plot F(x)
        plt.figure(figsize=(8, 5))
        plt.plot(x_vals, numerical_F_vals, label="CG $F$", linewidth=2)
        plt.xlabel("x")
        plt.ylabel(r"$F(x)$")
        plt.title(f"$F$ at t = {self.Tmax:.4f}")
        plt.legend()
        plt.grid(True)
        plt.tight_layout()
        plt.savefig(os.path.join(plot_dir, f"F_vs_exact_t_{self.Tmax:.4f}.png"))
        plt.close()

I obtain the following plots for Er:

and for F:

Clearly the plot for F is all wrong, and I can’t find a way to initialize F correctly.

I was wondering if I can get some help with this.

[connorjward]

Are you familiar with Firedrake's pre-existing matplotlib support: https://www.firedrakeproject.org/visualisation.html#plotting-with-matplotlib

Note that Function.at is deprecated practice and we recommend using a vertex-only-mesh instead: https://www.firedrakeproject.org/point-evaluation.html#primary-api-interpolation-onto-a-vertex-only-mesh

[afgalindo]

Thanks!!