Hello everyone! I am elaborating on the example provided at Elastic 3D beam structures — Computational Mechanics Numerical Tours with FEniCSx which implements a beam element model on a 1D mesh embedded in 3D using a mixed function space ((P^n, P^n), the first vector component for the domain’s displacement and the second vector component for the beam’s local rotation). Since my problem has no Dirichlet boundary conditions [my domain is kinda floating in outer space] I tried to make it solvable by setting a nullspace containing rigid translations and rotations. I took inspiration from the demo_elasticity.py example where we have
def build_nullspace(V: FunctionSpace):
"""Build PETSc nullspace for 3D elasticity"""
# Create vectors that will span the nullspace
bs = V.dofmap.index_map_bs
length0 = V.dofmap.index_map.size_local
basis = [la.vector(V.dofmap.index_map, bs=bs, dtype=dtype) for i in range(6)]
b = [b.array for b in basis]
# Get dof indices for each subspace (x, y and z dofs)
dofs = [V.sub(i).dofmap.list.flatten() for i in range(3)]
# Set the three translational rigid body modes
for i in range(3):
b[i][dofs[i]] = 1.0
# Set the three rotational rigid body modes
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.flatten()
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
b[3][dofs[0]] = -x1
b[3][dofs[1]] = x0
b[4][dofs[0]] = x2
b[4][dofs[2]] = -x0
b[5][dofs[2]] = x1
b[5][dofs[1]] = -x2
la.orthonormalize(basis)
basis_petsc = [
PETSc.Vec().createWithArray(x[: bs * length0], bsize=3, comm=V.mesh.comm) # type: ignore
for x in b
]
return PETSc.NullSpace().create(vectors=basis_petsc) # type: ignore
but in my case I have a mixed function space [displacements and local rotation] and my nullspace needs to include [I believe] translations and rotations for the displacement and translations for the rotation.
I modified the code above as follows
def _build_nullspace(W: fem.FunctionSpace) -> PETSc.NullSpace:
"""Create the operator nullspace.
This is used to impose average zero displacement and rotation to the solution.
"""
# Create vectors that will span the nullspace
block_size = W.dofmap.index_map_bs
local_size = W.dofmap.index_map.size_local
# AAA: We will need 9 vectors: six for the translations and rotations on the displacement,
# which is the first component of the space W, and three for the local rotation, which is the
# second component of the space W
nullspace_basis = [
la.vector(W.dofmap.index_map, bs=block_size, dtype=default_scalar_type)
for _ in range(9)
]
nullspace_basis_arrays = [b.array for b in nullspace_basis]
# Get dof indices for each subspace (x, y and z dofs)
# AAA: storing the dofs for the two subspaces separately
dofs_d = [W.sub(0).sub(i).dofmap.list.flatten() for i in range(3)]
dofs_r = [W.sub(1).sub(i).dofmap.list.flatten() for i in range(3)]
# Set the three translational rigid body modes
for i in range(3):
nullspace_basis_arrays[i][dofs_d[i]] = 1.0
# AAA: basis vectors 3, 4, 5 will be for the local rotation
for i in range(3, 6):
nullspace_basis_arrays[i][dofs_r[i-3]] = 1.0
# AAA: ok, now I have to set the rotations on the displacement.
# The main issue is that I cannot call W.tabulate_dofs_coordinates() so I
# have to collapse a subspace and do it there, to later map the coordinates back.
V, V_to_W = W.sub(0).collapse()
# AAA: I need to invert the map V_to_W in order to map the dofs back to the space W.
# For now I am creating a numpy array that maps the positions back so that I can use
# it for indexed/masked access conveniently. I put a value that is out of range in the missing
# entries so that if something goes wrong it will error out. Surely there's a better way to do this.
W_to_V = []
idx = 0
V_to_W_set = set(V_to_W)
for i in range(V_to_W[-1]+1):
if i in V_to_W_set:
W_to_V.append(idx)
idx += 1
else:
W_to_V.append(len(V_to_W)+1) # sketchy...
W_to_V = np.array(W_to_V)
# AAA: Now I can tabulate the dofs
x = V.tabulate_dof_coordinates()
# AAA: and get the dofmap, which refers to the space W
dofs_block = W.sub(0).dofmap.list.flatten()
# AAA: here I map the dofs back to the V space, dividing by three because
# dofs are grouped by three for a vector-valued function in 3D
idx_V = W_to_V[dofs_block] // 3
# AAA: Finally, I extract the dofs coordinates, picking only every third
# dof coordinate because they repeat three times
x0, x1, x2 = x[idx_V, 0][::3], x[idx_V, 1][::3], x[idx_V, 2][::3]
nullspace_basis_arrays[6][dofs_d[0]] = -x1
nullspace_basis_arrays[6][dofs_d[1]] = x0
nullspace_basis_arrays[7][dofs_d[0]] = x2
nullspace_basis_arrays[7][dofs_d[2]] = -x0
nullspace_basis_arrays[8][dofs_d[2]] = x1
nullspace_basis_arrays[8][dofs_d[1]] = -x2
la.orthonormalize(nullspace_basis)
basis_petsc = [
PETSc.Vec().createWithArray(
x[: block_size * local_size], bsize=3, comm=W.mesh.comm
)
for x in nullspace_basis_arrays
]
return PETSc.NullSpace().create(vectors=basis_petsc)
Now, incredibly enough, this horrible thing here seems to work, at least for a toy problem, at least in serial.
Here come my questions:
- Is my implementation actually correct?
- Is there a better way of doing this that doesn’t require sketchy manipulations of the dofs?
Thanks for any help!
Edit: Actually, it only works for P1-P1 elements so there’s definitely something off here.