Please see DOLFINx_MPC demos for the relevant version, As I believe this is due to you using the v0.8.x or main branch demos as a starting point: dolfinx_mpc/python/dolfinx_mpc/utils/mpc_utils.py at v0.7.2 · jorgensd/dolfinx_mpc · GitHub
i.e. using the nullspace function form v0.7.2 works fine:
def rigid_motions_nullspace(V: fem.FunctionSpaceBase):
"""
Function to build nullspace for 2D/3D elasticity.
Args:
V: The function space
"""
_x = fem.Function(V)
# Get geometric dim
gdim = V.mesh.geometry.dim
assert gdim == 2 or gdim == 3
# Set dimension of nullspace
dim = 3 if gdim == 2 else 6
# Create list of vectors for null space
nullspace_basis = [_x.vector.copy() for _ in range(dim)]
from contextlib import ExitStack
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]
dofs = [V.sub(i).dofmap.list.reshape(-1) for i in range(gdim)]
# Build translational null space basis
for i in range(gdim):
basis[i][dofs[i]] = 1.0
# Build rotational null space basis
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.reshape(-1)
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
if gdim == 2:
basis[2][dofs[0]] = -x1
basis[2][dofs[1]] = x0
elif gdim == 3:
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
la.orthonormalize(nullspace_basis)
assert la.is_orthonormal(
nullspace_basis, float(500 * np.finfo(_x.x.array.dtype).resolution)
)
return PETSc.NullSpace().create(vectors=nullspace_basis) # type: ignore
If I then add null_space.remove(uh.vector)
after solving the problem, I get an almost sensible solution.
The current issue is that you are mapping the bottom-most dof to itself.
What should the boundary condition be at that node?