Utilization of null space in elasticity problems

General
1

Hello everyone, I have utilized the null space while solving an elasticity problem in a 2D rectangular domain. I applied horizontal traction at the boundary where x=0. Below is the code：

from dolfin import *
import ufl
import numpy as np

# Create mesh and define function space
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

# Parameters
E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return 0.5*(grad(v) + grad(v).T)
def sigma(v):
    return  lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
class Leftboundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary
left = Leftboundary()

# Discrete traction
T_data = [0.01] * 11
T_coordinates = [[0, 0], [0, 0.1], [0, 0.2], [0, 0.3], [0, 0.4], [0, 0.5], [0, 0.6], [0, 0.7], [0, 0.8], [0, 0.9], [0, 1.0]]
T_space = FunctionSpace(mesh, "CG", 1)
T_load = Function(T_space)

vertex_mark = MeshFunction("size_t", mesh, 0, 0)
left.mark(vertex_mark, 1)
vertices = vertex_mark.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(T_space)

# Get all boundary coordinates
dof_coords = T_space.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array([dof_coords[boundary_dof] for boundary_dof in boundary_dofs])

points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(T_coordinates, 0)
distances = np.sum(np.square(points_A - points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    T_load.vector()[boundary_dofs[row]] = T_data[col]

boundary_mark = MeshFunction("size_t", mesh,1)
left.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

a = inner(sigma(u), eps(v)) * dx
Traction = ufl.as_vector((T_load, 0))
L = dot(Traction, v) * ds(1)
# Assemble system
A = assemble(a)
b = assemble(L)
# Solution Function
u = Function(V)
# Create Krylov solver
solver = PETScKrylovSolver("cg")
solver.set_operator(A)
# Create vector that spans the null space and normalize
null_vec = Vector(u.vector())
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")
# Create null space basis object and attach to PETSc matrix
null_space = VectorSpaceBasis([null_vec])
as_backend_type(A).set_nullspace(null_space)
null_space.orthogonalize(b)
solver.solve(u.vector(), b)

However, the runtime output is unclear, and I am unsure of where the issue lies:

Traceback (most recent call last):
  File "/home/dyfluid/Desktop/contact2DV/test.py", line 74, in <module>
    solver.solve(u.vector(), b)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

I hope to receive some assistance, thank you.

2

You are using the wrong nullspace, see:
https://bitbucket.org/fenics-project/dolfin/src/1c52e837eb54cc34627f90bde254be4aa8a2ae17/python/test/unit/la/test_nullspace.py?at=master#lines-26:60

4

Actually, I am not sure if null space can achieve this functionality: applying traction conditions to one edge of a 2D rectangular region and using null space to eliminate the translation and rotation of the rigid body. If it is possible, could you please make some modifications to my code below？ (this is my simple attempt, no errors occurred, but there are still significant displacements in the results). If it cannot be achieved, I will give up this method. Thank you.

from dolfin import *
import matplotlib.pyplot as plt

# Create mesh and define function space
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 2)

def build_nullspace(V, x):
    """Function to build nullspace for 2D/3D elasticity"""
    # 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.copy() for i in range(dim)]
    # Build translational null space basis
    for i in range(gdim):
        V.sub(i).dofmap().set(nullspace_basis[i], 1.0);
    # Build rotational null space basis
    if gdim == 2:
        V.sub(0).set_x(nullspace_basis[2], -1.0, 1);
        V.sub(1).set_x(nullspace_basis[2], 1.0, 0);
    elif gdim == 3:
        V.sub(0).set_x(nullspace_basis[3], -1.0, 1);
        V.sub(1).set_x(nullspace_basis[3],  1.0, 0);
        V.sub(0).set_x(nullspace_basis[4],  1.0, 2);
        V.sub(2).set_x(nullspace_basis[4], -1.0, 0);
        V.sub(2).set_x(nullspace_basis[5],  1.0, 1);
        V.sub(1).set_x(nullspace_basis[5], -1.0, 2);

    for x in nullspace_basis:
        x.apply("insert")
    # Create vector space basis and orthogonalize
    basis = VectorSpaceBasis(nullspace_basis)
    basis.orthonormalize()

    return basis

# Define boundaries
class Leftboundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary
left = Leftboundary()
boundary_mark = MeshFunction("size_t", mesh,1)
left.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

# Parameters
E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return 0.5*(grad(v) + grad(v).T)
def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(len(v))

# Define variational problem
T = Constant((0.1, 0))
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(T, u_) * ds(1)

A = assemble(a)
b = assemble(L)

# Create solution function
u = Function(V)
# Attach near nullspace to matrix
null_space = build_nullspace(V, u.vector())
as_backend_type(A).set_near_nullspace(null_space)
solve(A, u.vector(), b)
File("result/U.pvd") << u