Project function on custom vector space

dolfinx
2

I’m not really sure exactly what you’re trying to accomplish; however, it reads as though you’re aiming for some kind of Gram-Schmidt orthogonalisation process. Then you want to remove the rigid body mode components from a vector function. Perhaps the following helps, but please keep in mind I may have made mistakes:

from dolfin import *

mesh = UnitSquareMesh(32, 32)

V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)

# Some velocity function
u.interpolate(Expression(("-x[1]*x[1]", "x[0]*x[0]"), degree=4))

# Rigid body modes
rbms = [
    Constant((1.0, 0.0)),
    Constant((0.0, 1.0)),
    Expression(("-x[1]", "x[0]"), degree=1)
]
# List of functions to orthogonalise
v = list(interpolate(rbm, V) for rbm in rbms)

# GS Projection
def proj(u, v):
    res = assemble(inner(u, v)*dx)/assemble(inner(u, u)*dx) * u.vector()
    return res

# GS orthogonalisation
def ortho(v):
    xi = [None]*len(v)
    xi[0] = v[0]
    for j in range(1, len(xi)):
        xi[j] = Function(V, v[j].vector() - sum(proj(xi[i], v[j]) for i in range(j)))
    return xi

# Orthogonalised vectors
xi = ortho(v)

# Orthonormalised vector basis
e = [Function(V, xi_.vector()/assemble(inner(xi_, xi_)*dx)**0.5) for xi_ in xi]

print("orthonormalisation test:")
for i in range(len(xi)):
    for j in range(i+1):
        print(f"inner(e[{i}], e[{j}])*dx {assemble(inner(e[i], e[j])*dx)}")

u_star = u.vector() - sum(proj(e_, u) for e_ in e)
u_star = Function(V, u_star)

print(f"u norm {u.vector().norm('l2')}, u_star norm {u_star.vector().norm('l2')}")
print(f"orthogonalisation of u_star with rigid body modes test:")
for j in range(len(rbms)):
    print(f"(rbms[{j}], u_star) = {assemble(inner(u_star, rbms[j])*dx)}")

which outputs the following for me:

orthonormalisation test:
inner(e[0], e[0])*dx 1.0
inner(e[1], e[0])*dx 0.0
inner(e[1], e[1])*dx 1.0
inner(e[2], e[0])*dx 1.1926223897340549e-18
inner(e[2], e[1])*dx -3.2526065174565133e-19
inner(e[2], e[2])*dx 1.0
u norm 21.35620361537434, u_star norm 3.6983169052517852
orthogonalisation of u_star with rigid body modes test:
(rbms[0], u_star) = -1.4907779871675686e-19
(rbms[1], u_star) = 1.0394788328704774e-17
(rbms[2], u_star) = -9.378348791999613e-18