Orthonormal basis of FunctionSpace

Linear Algebra
1

Hi,

I base myself in the context of generalized polynomial chaos expansions for uncertainty quantification. To this end, I can prove that the coefficients of the solutions of certain elliptic PDEs in an orthonormal basis of the function space are not changing too much (this statement can be made more precised, but that’s not the topic here).

In order to have nice numerical experiments with this, I wanted to implement this. Here is the problem I am facing now: classical finite element basis are not orthonormal. I found the VectorSpaceBasis in the Dolfin.la class which has a Gram-Schmidt function, but this only works – as far as I understand – on discrete vectors. Is there something similar directly useable in Fenics?

from fenics import *
import numpy as np

# First create a 2D mesh (better for understanding)
nb_points = 10
msh = UnitSquareMesh(nb_points,nb_points)

# Function spaces that will be useful to us
V = FunctionSpace(msh, 'Lagrange', 1)

# Define boundary conditions
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
w = TrialFunction(V)
v = TestFunction(V)

# Some parameters for playing
Aparam = 1.01 # Slightly perturbed problem
Lparam = -6.0
# And here comes the problem
x = SpatialCoordinate(msh)
A = Constant(Aparam) * inner(nabla_grad(w), nabla_grad(v)) * dx
L = Constant(Lparam) * v * dx

# Create goal-functional for error estimation
u = Function(V)
phi1 = Function(V)

nb_dim = V.dim()
FEbasis = []

for i in range(nb_dim): 
    e1 = np.zeros(nb_dim)
    e1[i] = 1
    phi1.vector().set_local(e1)
    FEbasis.append(phi1)

local_basis = FEniCS.Basis(FEbasis).orthonormalize() # Does something like this exists? 

solve(A==L, u, bc)
# Analyze the coefficients of the slightly perturbed solution w.r.t ONB
inner_prod[0] = assemble(inner(u, local_basis[0])*dx)
inner_prod[1] = assemble(inner(u, local_basis[1])*dx)
# etc ...

I would rather not have to compute the whole Gram-Schmidt process by myself (although the time I have spent trying to find ways around would have probably been sufficient to implement it … )

Thanks in advance for any help!

2

Hi, is it the L^2 inner product that you wish to use for othogonalization of the basis?

3

Thanks for this reply, I should I have mentioned it indeed.

Yes, an L^2 orthonormalization will be perfect as a first thing. Later I might have to adapt a wee bit, but likely in a L^2 + L^2(\nabla), but this shouldn’t be a major issue afterwards.