# Verifying Null Space in Mixed Space Linear Elasticity Problems

Hi all,
I am working on a simple elasticity problem and I would like to add a new variable to that which should be equal to the displacment gradient:
here is the variational formulation:

πππ: the stress tensor, πππ : the strain tensor, and π’π : the displacement field. Ξ¨jk is the new variable that should equals the π’π,π (du/dxj). I used Lagrange multiplier to force this relation.

Since I am using itertaive solvers, I need to set up the null space for this variational form.
I think I have to provide 18 vectors that span the null space:( 6 for rigid body motions, 6 for gradient of the displacement field π’π,π,and 6 arbitrary functions for πΌππ).

I followed the post here to set up the null space and was able to do that. But I donβt know how to make sure this is the correct null space for my problem. What are the option to make sure that I have set up the correct nulll space?
I tried to test it with

``````assert nullspace.test(A)
``````

but I got the assertaion error.

Here is MWE:

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx.fem import (Function, functionspace, dirichletbc, form, locate_dofs_topological, Constant)
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.mesh import locate_entities_boundary, create_box, CellType, create_unit_square
from ufl import (Identity, TestFunctions, TrialFunctions, dx, indices, as_tensor)
from basix.ufl import element, mixed_element
from dolfinx import la
from dolfinx import default_scalar_type
import dolfinx
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
F_star = Fs
# ........................................ Parameters
E = 1 * 152e9 * C_star
nu = 0.33
c1 = E * nu / ((1 + nu) * (1 - 2 * nu))
c2 = E / (2 * (1 + nu))
h = 0.76e-3 / Fl
H = h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],[1, 1, 1], cell_type=CellType.hexahedron)
# ......................................... Define Elements
v = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
s = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Displacement gradient Element
lm = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Lagrange multipliers Element
# ......................................... Define mixed element
MIXEL = mixed_element([v, s, lm])
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
# ....................................... Mark boundary regions
tol = 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)
# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
clamp = Function(S0)
clamp.x.array[:] = 0
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
# .........Fixed BC
clamp_dofs = locate_dofs_topological((Space.sub(0), S0), mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(clamp, clamp_dofs, Space.sub(0))
bc = [bc0]
# ......................................... Define tensor
delta = Identity(len(u))
i, j, k = indices(3)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

def StressT(u):
return as_tensor((c1 * StrainT(u)[i, i] * delta[j, k] + 2 * c2 * StrainT(u)[j, k]), [j, k])
return as_tensor(gradu[j, k] - u[k].dx(j), (j, k))
# ......................................... Define Variational Formulation
a = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a += Constraint(u, gradu)[j, k] * del_lamu[j, k] * dx
a += Constraint(del_u, del_gradu)[j, k] * lamu[j, k] * dx
#.................................................Assembling
bilinear_form = form(a)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
# ################################### Apply Nullspace
def build_nullspace(V: Space):
bs = Space.dofmap.index_map_bs  # number of block
length0 = Space.dofmap.index_map.size_local  # number of dof in the cell
basis = [la.vector(Space.dofmap.index_map, bs=bs, dtype=PETSc.ScalarType) for i in range(18)]  # 18 is the number of basis we need
b = [b.array for b in basis]
# Get dof indices for each subspace (x, y and z dofs)
dof_u = [Space.sub(0).sub(i).dofmap.list.flatten() for i in
range(3)]  # creates dof associated with each nodal variables
# Set the three translational rigid body modes
for i in range(3):
b[i][dof_u[i]] = 1.0
# Set the three rotational rigid body modes
S0, S0_to_Space = Space.sub(0).collapse()
x = S0.tabulate_dof_coordinates()
dofs_block = S0.dofmap.list.flatten()
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
b[3][dof_u[0]] = -x1
b[3][dof_u[1]] = x0
b[4][dof_u[0]] = x2
b[4][dof_u[2]] = -x0
b[5][dof_u[2]] = x1
b[5][dof_u[1]] = -x2
##################################
dof_gu = [Space.sub(1).sub(i).dofmap.list.flatten() for i in range(9)]
b[6][dof_gu[0]] = 1.0
b[7][dof_gu[4]] = 1.0
b[8][dof_gu[8]] = 1.0
b[9][dof_gu[5]] = 1.0
b[9][dof_gu[7]] = -1.0
b[10][dof_gu[2]] = 1.0
b[10][dof_gu[6]] = -1.0
b[11][dof_gu[1]] = 1.0
b[11][dof_gu[3]] = -1.0
#######################
dof_lam = [Space.sub(2).sub(i).dofmap.list.flatten() for i in
range(9)]  # creates dof associated with each nodal variables

b[12][dof_lam[0]] = 1.0
b[13][dof_lam[1]] = 1.0
b[13][dof_lam[3]] = 1.0
b[14][dof_lam[2]] = 1.0
b[14][dof_lam[6]] = 1.0
b[15][dof_lam[4]] = 1.0
b[16][dof_lam[5]] = 1.0
b[16][dof_lam[7]] = 1.0
b[17][dof_lam[8]] = 1.0

_basis = [x._cpp_object for x in basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_basis)

basis_petsc = [
PETSc.Vec().createWithArray(x[: bs * length0], bsize=1, comm=Space.mesh.comm)  # type: ignore
for x in b
]
return PETSc.NullSpace().create(vectors=basis_petsc)

nullspace = build_nullspace(Space)
assert nullspace.test(A)
``````

β

β

1 Like

I would suggest starting with a slightly simpler problem. Start with checking that you can build the nullspace for a simple elastic problem, then ensure if you make a mixed formulation of two elasticity problems (they can be equivalent) that the nullspace creating logic is still correct.

Itβs quite alot of work for anyone to first ensure that you have implemented a sensible variational form, and then in turn as them to debug 18 null-spaces for you. I would start with ensure that rigid body motions are well defined before going further.

Thank you so much for your answer. I am quite sure the variational formulation is correct, as I have tested it with a direct solver. I believe the issue lies with the null space.
I really appreciate your idea of using a mixed formulation for two elasticity problems. However, could you please explain what you mean by that? Do you mean defining two different vector elements and then mixing them? If so, there wouldnβt be any coupling between them. Is that correct?

As you suggested, as the first step, I attempted this for a simple elasticity problem. I considered a cube with traction applied to both the top and bottom surfaces (pure Neumann boundary conditions):I successfully set up the null space without any errors, and everything appeared to be functioning correctly. However, I observed that there was no difference in the results with or without using the null space. I mean, it seems there was no need to build the nullspace for this problem. So:

• How can I ensure that I have correctly built the null space? What should I check?

Here is the code :

``````import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological, Expression)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
from dolfinx import la
#
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# # # ....................
C_star = Fl ** 2 / Fs
# F_star = Fs
# -------------------------Material Properties
E = 1 * 130e9 * C_star
nu = 0.28
laa = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
h = 0.76e-3 / Fl
H= h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],[15, 15, 15], cell_type=CellType.hexahedron)
Force = 200.00000 / Fs
Area1 = h**2 / Fl ** 2  # Bottom
Area2 = h**2 / Fl ** 2
tr_S1 = Constant(mesh, ScalarType((0, 0, Force / Area1)))  # Bottom
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
# ......................................... Define mixed element
MIXEL = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# # ...................................... BCs
bc = []
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l  = indices(4)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l], (i, j, k, l))  # Elasticity tensor
def StressT(u):
return as_tensor((Elas(laa, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
# ......................................... Define Variational Formulation
a = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
L = tr_S2[i] * del_u[i] * ds(2) + tr_S1[i] * del_u[i] * ds(1)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
bn = assemble_vector(linear_form)
apply_lifting(bn, [bilinear_form], [bc])
set_bc(bn, bc)
##################
dtype = PETSc.ScalarType  # type: ignore
def build_nullspace(V: Space):
"""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

_basis = [x._cpp_object for x in basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_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

ns = build_nullspace(Space)
assert ns.test(A)
A.setNearNullSpace(ns)
################
Result = Function(Space)
def monitor(ksp, it, rnorm):
print(f"iteration: {it} rnorm: {rnorm:.3e}")
# Set solver options
solver = PETSc.KSP().create(mesh.comm)  # type: ignore
opts = PETSc.Options()  # type: ignore
opts["ksp_type"] = "gmres"
opts["ksp_rtol"] = 1.0e-15
opts["ksp_atol"] = 1.0e-15
opts["pc_type"] = "gamg"
solver.setFromOptions()
solver.setOperators(A)
print(solver.getConvergedReason())
solver.setMonitor(monitor)
solver.solve(bn, Result.vector)
print(solver.getConvergedReason())
Result.x.scatter_forward()
##################
V_u = functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
disp = Function(V_u)
disp_expr = Expression(Result * Fl, V_u.element.interpolation_points())
disp.interpolate(disp_expr)

with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [disp]) as vtx:
vtx.write(0.0)
``````

Nullspaces are tested as you have already illustrated:

If this test fails it means that the built nullspaces is not a null-space for the assembled operator `A`.

Yes, that would be the initial idea, to ensure that you have understood how to build up null-spaces for a mixed element.