Error: Unable to orthonormalize vector basis. Reason: Vector space has linear dependency. Where: This error was encountered inside Vector Space Basis.cpp

General
10

Because the tutorial is very out of date.
Local project shouldn’t consume alot of memory, as it only assembles local matrices, and does local inversions, Which should be quicker and less expensive than a normal projection.

I cannot reproduce any of your errors, as you haven’t supllied the msh file that you are using. Here is a version with the original mesh you posted:

from dolfin import *
from mshr import *

def build_nullspace(V, x):

    gdim = V.mesh().geometry().dim()
    assert gdim == 2 or gdim == 3

    dim = 3 if gdim == 2 else 6

    nullspace_basis = [x.copy() for i in range(dim)]

    for i in range(gdim):
        V.sub(i).dofmap().set(nullspace_basis[i], 1.0);

    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")

    return VectorSpaceBasis(nullspace_basis)

cylinder = Cylinder(Point(0.5, 0.5, 0), Point(0.5, 0.5, 2), 0.25, 0.25)
cube = Box(Point(0, 0, 0),Point(1, 1, 1))
domain = cube - cylinder
mesh = generate_mesh(domain, 40)


E = Constant(200e3) #E = 200 GPa 
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return lmbda*tr(epsilon(u))*Identity(3) + 2*mu*epsilon(u)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

f = Constant((0., 0., 0.))
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), epsilon(v))*dx
L = inner(f, v)*dx

class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0., 0.1)
class right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1., 0.1)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

left().mark(boundaries, 1)
right().mark(boundaries, 2)

bc_left = DirichletBC(V.sub(0), Constant(1.) , boundaries, 1)
bc_right = DirichletBC(V , Constant((0., 0., 0.)), boundaries, 2)

bcs = [bc_left, bc_right]

A, b = assemble_system(a, L, bcs)

u = Function(V)

null_space = build_nullspace(V, u.vector())

as_backend_type(A).set_near_nullspace(null_space)

pc = PETScPreconditioner("petsc_amg")

PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")

PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)

solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = True

solver.set_operator(A);

solver.solve(u.vector(), b);


def local_project(f, V):
    u = TrialFunction(V)
    v = TestFunction(V)
    a_proj = inner(u, v)*dx
    b_proj = inner(f, v)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    u = Function(V)
    solver.solve_local_rhs(u)
    return u

W = TensorFunctionSpace(mesh, "Discontinuous Lagrange", 0)
stress = local_project(sigma(u), W)
with XDMFFile("stress.xdmf") as xdmf:
    xdmf.write(stress)