Project function on custom vector space

dolfinx
1

Hello,

Let’s say I have a function u that represents a 3D CG1 displacement field. What I would like to do is to extract the rigid body component of u, e.g. break down u = u_0 + u_rigid.

I know how to construct the nullspace for 3D elasticity :

def nullspace_elasticity_3D(V):
    """Function to build null space for 3D elasticity"""

    # Create list of vectors for null space
    index_map = V.dofmap.index_map
    nullspace_basis = [dolfinx.cpp.la.create_vector(index_map, V.dofmap.index_map_bs) for i in range(6)]

    with ExitStack() as stack:
        vec_local = [stack.enter_context(x.localForm()) for x in nullspace_basis]
        basis = [np.asarray(x) for x in vec_local]

        # Dof indices for each subspace (x, y and z dofs)
        dofs = [V.sub(i).dofmap.list.array for i in range(3)]

        # Build translational null space basis
        for i in range(3):
            basis[i][dofs[i]] = 1.0

        # Build rotational null space basis
        x = V.tabulate_dof_coordinates()
        dofs_block = V.dofmap.list.array
        x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
        basis[3][dofs[0]] = -x1
        basis[3][dofs[1]] = x0
        basis[4][dofs[0]] = x2
        basis[4][dofs[2]] = -x0
        basis[5][dofs[2]] = x1
        basis[5][dofs[1]] = -x2

    # Create vector space basis and orthogonalize
    basis = dolfinx.cpp.la.VectorSpaceBasis(nullspace_basis)
    basis.orthonormalize()

    _x = [basis[i] for i in range(6)]
    nsp = PETSc.NullSpace().create(vectors=_x)
    return nsp

So what I wanted to do, is take the VectorSpaceBasis defined in this function, create a VectorFunctionSpace from it, and then project u on this space to get what I want.

However I don’t know if it is possible to create a VectorFunctionSpace from a VectorSpaceBasis ?
If not what would be the prefered way of doing what I want ?

Thanks!

1 Like
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
3

Thanks nate, that seems to do the trick. For some context, I was trying to compare results of FEM simulation with measured data. The measured displacement field was measured from a video and contains a rigid motion. However my simulation is expressed in terms of a pure Neumann elasticity PDE, so I constrained my solution to have no rigid motion. Hence the need for me to remove the rigid motion from the measured displacement.

4

Sorry for reviving an older thread, but I just wanted to share the dolfinx version of this code in case anyone finds it useful. Thanks for this piece of code in the first place @nate. Hopefully others find this useful if they need to do any orthogonalization w.r.t. a custom vector space.

from dolfinx import fem,mesh
from mpi4py import MPI
from ufl import inner,Measure,as_vector,SpatialCoordinate
from petsc4py import PETSc

nx=32
ny=32
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)

V = fem.VectorFunctionSpace(domain, ("CG", 1))
u = fem.Function(V)
x = SpatialCoordinate(domain)

#Some random expression
u_expr = fem.Expression(as_vector([-x[1]*x[1],x[0]*x[0]]),V.element.interpolation_points())
u.interpolate(u_expr)

dx = Measure("dx",domain=domain)

#Rigid Body Modes expression (2D)
rbms = [
    fem.Expression(fem.Constant(domain,PETSc.ScalarType((1.0,0.0))),V.element.interpolation_points()),
    fem.Expression(fem.Constant(domain,PETSc.ScalarType((0.0,1.0))),V.element.interpolation_points()),
    fem.Expression(as_vector([-x[1],x[0]]),V.element.interpolation_points()),
]

# List of functions to orthogonalise
vx = fem.Function(V)
vy = fem.Function(V)
vr = fem.Function(V)
vx.interpolate(rbms[0])
vy.interpolate(rbms[1])
vr.interpolate(rbms[2])

v = list((vx,vy,vr))

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

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

# Orthogonalised vectors
xi = ortho(v)

# Orthonormalised vector basis
e = [fem.Function(V) for i in range(len(v))]
for i,xi_ in enumerate(xi):
    e[i].vector.array = xi_.vector.array/fem.assemble_scalar(fem.form(inner(xi_, xi_)*dx))**0.5

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

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

print(f"u norm {u.vector.norm(2)}, u_star norm {u_star.vector.norm(2)}")
print(f"orthogonalisation of u_star with rigid body modes test:")
for j in range(len(v)):
    print(f"(rbms[{j}], u_star) = {fem.assemble_scalar(fem.form(inner(u_star, v[j])*dx))}")
2 Likes
5

Hi @jkrokowski and @nate , thanks for this bit of code.

I have a question about the orthogonalisation:

Could you explain in the proj function why it is necessary to assemble the form fem.assemble_scalar(fem.form(inner(u, v)*dx)), rather than something like np.dot(u.x.array,v.x.array)?

6

The difference here is between the l2 and L2 norm of two functions from a discrete function space.

the numpy dot product does not take into account that the mesh size might be different between the different degrees of freedom.

1 Like
7

Thanks Dokken, makes perfect sense.

So, the l2 norm might go wrong if the distribution of nodes in a mesh is not completely regular, whereas L2 will always be correct?